A tiny distributed genetic programming engine (work in progress).
A few years ago I read Genetic Programming - On The Programming of
Computers by Means of Natural Selection by John Koza, published in 1992 by MIT.
I wasted many hours building GP systems in Python and C++ before I really
understood why he used Lisp.
Now I'm learning Scheme and functional programming, and I decided to write up
my notes and my code as I build a distributed GP system. Just for fun.
Let's say we have the following data points, and we want to grow a reasonable definition of the function f:
x f(x)
------+-------
-15 | 862
1 | 14
10 | 437
42 | 7189
100 | 40307
These data points happen to be for the function
f(x) = 4x² + 3x + 7, but we want to 'grow' a suitable
function, without using this knowledge. The first step in setting up a GP
experiment is to decide on a reasonable set of functions and terminals to use
in the genetic programs. A good base set of functions for problems like this
is the four arithmetic operators (although we use a special 'protected divide'
which allows division by zero, returning zero in this case, traditionally
indicated by the symbol %, from the Koza book),
but could also include trigonometric
and logarithm functions. When specifying the functions, we
need to state the number of parameters they take. For terminals, we need
to include the variable
x, and then the special symbol random-int which will
produce arbitrary numbers from 0 to 10 (Koza calls this the 'ephemeral random
constant'). Let's set up some parameters to play with by typing straight into
the REPL (everything shown in red and green is just throw-away test interaction
with the Scheme interpreter - whereas the blue and orange are extracts of
source code, which can of course also be entered straight into the REPL but
in real use would be loaded from a library):
> (define my-functions '#((+ . 2) (- . 2) (% . 2) (* . 2)))> (define my-terminals '#(random-int x))> (define my-min-depth 2)> (define my-max-depth 4)
At the heart of the GP system is the individual, which consists of a
program tree, and a fitness value. These are implemented using a Scheme pair,
so cons, car and cdr would be enough
to work with them, but it's clearer with explicit functions:
; Creates an individual from a program and a fitnes.
(define (make-individual program fitness)
(cons program fitness))
; Returns the program tree of an individual.
(define (individual-program individual)
(car individual))
; Returns the fitness of an individual.
(define (individual-fitness individual)
(cdr individual))
The first thing we want to do is create new random programs for the initial population:
; Builds a random program from the provided vectors of functions and
; terminals.
(define (make-random-program functions terminals min-depth max-depth)
(define (iter depth)
(if (and (< depth max-depth)
(or (< depth min-depth) (zero? (random:random-integer 2))))
; generate a function application
(let* ([pair (vector-ref functions
(random:random-integer (vector-length
functions)))]
[function (car pair)]
[num-args (cdr pair)]
[arguments (let loop ([n num-args])
(if (zero? n)
'()
(cons (iter (+ 1 depth)) (loop (- n 1)))))])
(cons function arguments))
; generate a terminal node
(let ([terminal (vector-ref
terminals
(random:random-integer (vector-length terminals)))])
(cond [(eq? terminal 'random-real)
(- (* 10 (random:random-real)) 5.0)]
[(eq? terminal 'random-int)
(random:random-integer 10)]
[else terminal]))))
(iter 0))
We can test that this works in the REPL:
> (make-random-program my-functions my-terminals my-min-depth my-max-depth)
(+ (* x (+ (- x 9) (+ x 7))) (* 2 x))
Next we need a way to produce a whole population of them:
; Creates a list of random (program . fitness) pairs from the list of
; functions and terminals provided.
(define (make-population size functions terminals min-depth max-depth)
(define (iter n output)
(if (zero? n)
output
(iter (- n 1)
(cons (make-individual (make-random-program
functions terminals
min-depth max-depth)
#f)
output))))
(iter size '()))
We can show this working in the REPL:
> (make-population 10 my-functions my-terminals my-min-depth my-max-depth)
(((* (% (% (% x x) (% x 7)) (+ 2 x)) (+ 5 9))
.
#f)
((* (% x (+ x x)) (* (- 2 (* 9 x)) (* (% 8 x) 3)))
.
#f)
((% (* (+ (- 1 x) 4) (- (+ x x) (- x 3))) (- x 9))
.
#f)
((- (* 0 (+ 6 (+ x 0))) (- (- (% 6 7) 5) 4))
.
#f)
((* (- (* (* 4 9) (% x 3)) x) (% x 3)) . #f)
((% (- 2 (* (* 2 7) (- x 9))) (+ 3 x)) . #f)
((+ (* (% (% 0 x) (+ 5 3)) (+ (+ x 8) (+ x 5)))
(% x (+ 8 6)))
.
#f)
((- (+ 9 (+ (+ 6 x) x)) (+ (- x 8) (% (% x x) x)))
.
#f)
((- (* 4 (% x (% x 3))) (- (+ x 0) 8)) . #f)
((+ (* x (+ (- x 9) (+ x 7))) (* 2 x)) . #f))
Those #f values in each pair are the fitnesses, which have not
yet been decided. Next, we need a way to evaluate all of these programs!
; Computes and sets the fitness of each individual, and sorts into fitness
; order (fittest at the beginning).
(define (evaluate-population population fitness-function)
(define (iter input output)
(if (null? input)
output
(let ([program (caar input)]
[fitness (cdar input)])
; we only compute the fitness if it was not already known (#f)
(if (eq? fitness #f)
(iter (cdr input)
(cons (make-individual program (fitness-function program))
output))
(iter (cdr input)
(cons (make-individual program fitness)
output))))))
(list-sort (lambda (left right) (> (cdr left) (cdr right)))
(iter population '())))
To test this we'll need a fitness function, which is not part of the GP engine itself but part of the configuration for the experiment, so I'm going to paste it into the REPL:
> (define (my-fitness-test individual)
(define function (eval `(lambda (x) ,individual) (interaction-environment)))
(/ 1.0 (+ 1 (abs (- (function -15) 862))
(abs (- (function 1) 14))
(abs (- (function 10) 437))
(abs (- (function 42) 7189))
(abs (- (function 100) 40307)))))
The function returns 1 / (1 + error) where error is the sum of the how far wrong the candidate function was at each of the sample points. So
if we hit exactly the right function, the result would be 1.0, and if we were
very wrong, we'd have a number close to 0.0.
Notice that the fitness function is responsible for taking the program
tree and converting it into a Scheme function, which it does wrapping it in
lambda so that the name x is a function
parameter. This can't be done by Revolver because it doesn't know which of
the terminals should be function parameters, and in which order. Now let's
create a population, take a look at it, and evaluate it:
> (define my-generation-0 (make-population 10 '#((+ . 2) (- . 2) (% . 2) (* . 2)) '#(2 3 x) 2 4))
> my-generation-0
(((* (* (- (* 3 x) 3) (% (% 2 x) 3)) (- (% x 2) (- 3 3))) . #f)
((- (% x (% (* 3 2) (- 2 x))) (% (* 2 (- x 3)) 3)) . #f)
((+ (* 2 (% (% x 2) (% 3 2))) (+ (- 2 3) (+ 3 3))) . #f)
((+ (+ 2 (+ (- 2 2) 2)) (+ (+ 3 (- x 3)) 3)) . #f)
((- (+ (* (- 2 2) (* 2 3)) (% (+ x 3) (- 3 2))) (- (- (% x x) 2) (+ (+ x 2) (+ x x)))) . #f)
((+ (+ 2 (+ 2 (- 3 x))) (+ (+ (- 3 x) 3) 2)) . #f)
((+ (- 3 x) (+ 3 x)) . #f)
((* (- x 2) (- (- (% 2 x) x) x)) . #f)
((+ (+ (* x (% 3 2)) 2) (- 2 3)) . #f)
((% (* (- 2 (% 3 x)) (% x x)) (- 3 x)) . #f))
> (evaluate-population my-generation-0 my-fitness-test)
(((- (+ (* (- 2 2) (* 2 3)) (% (+ x 3) (- 3 2))) (- (- (% x x) 2) (+ (+ x 2) (+ x x)))) . 2.073484282989135e-05)
((+ (+ (* x (% 3 2)) 2) (- 2 3)) . 2.0576978476480515e-05)
((+ (+ 2 (+ (- 2 2) 2)) (+ (+ 3 (- x 3)) 3)) . 2.0560478647942926e-05)
((* (* (- (* 3 x) 3) (% (% 2 x) 3)) (- (% x 2) (- 3 3))) . 2.0543583211783798e-05)
((+ (* 2 (% (% x 2) (% 3 2))) (+ (- 2 3) (+ 3 3))) . 2.0536832809644097e-05)
((+ (- 3 x) (+ 3 x)) . 2.050020500205002e-05)
((% (* (- 2 (% 3 x)) (% x x)) (- 3 x)) . 2.048731521543056e-05)
((+ (+ 2 (+ 2 (- 3 x))) (+ (+ (- 3 x) 3) 2)) . 2.0403582869151823e-05)
((- (% x (% (* 3 2) (- 2 x))) (% (* 2 (- x 3)) 3)) . 1.9661430172430744e-05)
((* (- x 2) (- (- (% 2 x) x) x)) . 1.3806001382020183e-05))
Now we have the mechanisms to create the first generation and evaluate its fitness, so we now need a way to produce the next generation. The basic mechanisms for this are cross-breeding and duplication of programs. For cross-breeding we'll need a couple of helper functions for measuring and building program trees:
; Counts the atoms in a tree of lists and atoms.
(define (tree-atoms tree)
(cond [(null? tree) 0]
[(pair? tree) (+ (tree-atoms (car tree)) (tree-atoms (cdr tree)))]
[else 1]))
> (tree-atoms 'x)
1
> (tree-atoms '(+ 1 2)
3
> (tree-atoms '(+ x (* y z))
5
; Retrieves a node from a tree by index number, with a numbering scheme
; as shown:
;
; (+ x (* y z)) -> + -> 0
; / \ / \
; x * 1 2
; / \ / \
; y z 3 4
(define (tree-ref tree index)
; TODO find out how to do this without using a mutable counter to number
; the nodes (I have found only a horrible inefficient way)
(define count 0)
(call/cc
(lambda (return)
(define (parameters node)
(cond [(null? node) #f]
[else (test (car node))
(parameters (cdr node))]))
(define (test node)
(if (= count index) (return node))
(set! count (+ 1 count))
(if (pair? node) (parameters (cdr node))))
(test tree)
(error "TREE-REF -- index out of range"))))
> (tree-ref '(+ x (* y z)) 0)
(+ x (* y z))
> (tree-ref '(+ x (* y z)) 1)
x
> (tree-ref '(+ x (* y z)) 2)
(* y z)
> (tree-ref '(+ x (* y z)) 3)
y
> (tree-ref '(+ x (* y z)) 4)
z
> (tree-ref '(+ x (* y z)) 5)
ERROR: TREE-REF -- index out of range
; Copies a tree, inserting 'guest' into the place identified by 'index'.
;
; + +
; / \ % / \
; x <*> + / \ = x %
; / \ x 42 / \
; y x x 42
;
; TODO as above, I need to find out how to do this without the mutable count
(define (tree-combine host index guest)
(define count 0)
(define (parameters node)
(if (null? node) '()
(cons (copy (car node))
(parameters (cdr node)))))
(define (copy node)
(cond [(= count index)
(set! count (+ 1 count))
guest]
[else
(set! count (+ 1 count))
(if (pair? node)
(cons (car node) (parameters (cdr node)))
node)]))
(copy host))
> (tree-combine '(+ x (* y z)) 2 '(% x 42))
(+ x (% x 42))
Cross-breeding two individual programs is then as simple as:
; Creates random offspring by combining two program trees.
(define (tree-random-recombine male female)
(let* ([male-index (random:random-integer (tree-atoms male))]
[female-index (random:random-integer (tree-atoms female))]
[male-fragment (tree-ref male male-index)])
(tree-combine female female-index male-fragment)))
The next problem is how to select individuals for cross breeding. Koza proposed various methods but the most straightforward is to choose individuals with probability fitness divided sum of fitness. So if there were three individuals a, b, c and which had fitness 0.8, 0.6 and 0.6 respectively, then a should have a 0.8 / 2.0 or 40% chance of being selected, and b and c should both have a 30% chance. To do that efficiently we first build an index of fitness which we can use to pick individuals:
; Builds a vector of cumulative sums. This is useful for throwing a dart
; to pick an individual with fitness probability.
(define (make-fitness-index population)
(define (iter input output sum)
(cond [(null? input) output]
[else (iter (cdr input)
(cons sum output)
(+ sum (individual-fitness (car input))))]))
(list->vector (reverse (iter population '() 0.0))))
> (make-fitness-index '((a . 0.8) (b . 0.6) (c . 0.6)))
#(0 0.8 1.4)
This index provides us with a fast way to select an individual. We simply pick a random number between 0 and the sum of all fitness, and then we find position of the right-most value that is lower than our number. The position identifies the individual by index number. We therefore need a binary search:
; Finds a record which has a value equal to or immediately
; preceding a given value.
; get-ref - a function that can get a record by index
; size - the size of the set of records
; value - the search value
; precedes? - a function which can compare keys
(define (search get-ref size value precedes?)
(let loop ([start 0]
[stop (- size 1)])
(if (< stop start)
(if (< stop 0) #f stop)
(let* ([mid-point (quotient (+ start stop) 2)]
[key (get-ref mid-point)])
(cond [(precedes? value key) (loop start (- mid-point 1))]
[(precedes? key value) (loop (+ mid-point 1) stop)]
[else mid-point])))))
> (define my-vec #(0 0.8 1.4))
> (define (my-get n) (vector-ref my-vec n))
> (search my-get 3 0.3 <)
0
> (search my-get 3 0.8 <)
1
> (search my-get 3 1.5 <)
2
We finally have enough machinery to evolve a new generation:
; Breeds a new population. The new population is created by cloning and
; crossing invididuals in the source population. The population must be
; evaluated and sorted (fittest first).
(define (breed-population population crossover-probability)
(let* ([fitness-index (make-fitness-index population)]
[individuals (vector-length fitness-index)]
[population-vector (list->vector population)]
[crossovers (floor (* individuals crossover-probability))]
[copies (- individuals crossovers 1)]
[get-sum (lambda (n) (vector-ref fitness-index n))]
[total-fitness (vector-ref fitness-index (- individuals 1))])
; Selects an individual with fitness probability.
(define (select-individual)
(let* ([dart (* (random:random-real) total-fitness)]
[index (search get-sum individuals dart <)])
(vector-ref population-vector index)))
; Builds a list of 'count' combined individuals. Their fitness is set
; to #f because it will need to be recomputed.
(define (combine output count)
(if (zero? count) output
(let* ([male (individual-program (select-individual))]
[female (individual-program (select-individual))]
[offspring (tree-random-recombine male female)])
(combine (cons (make-individual offspring #f) output)
(- count 1)))))
; Builds a list of 'count' copies individuals (preserving their fitness
; so it will not be recomputed in the next generation).
(define (copy output count)
(if (zero? count) output
(copy (cons (select-individual) output)
(- count 1))))
; We want a list of combined and copied individuals, with the fittest
; individual thrown in at the front for good measure.
;
(cons (car population)
(append (copy '() copies)
(combine '() crossovers)))))
All we need now is the iteration machinery to evolve successive generations. Here's a simple version that does it and gives you the results:
; Evolves a population for 'generations' generations, or until the target
; fitness is reached, and returns the result.
(define (evolve population fitness-function generations crossover-probability
target-fitness)
(let loop ([generations generations]
[population population])
(define evaluated (evaluate-population population fitness-function))
(show-best evaluated)
(cond [(zero? generations) evaluated]
[(<= target-fitness (individual-fitness (car evaluated)))
evaluated]
[else
(loop (- generations 1)
(breed-population evaluated crossover-probability))])))
This code was mostly developed on DrScheme, which provides a nice testing and debugging environment for Scheme newbies like me. But for production GP runs, we need the fastest Scheme we can get our hands on. The word on the street is that Ikarus is very fast. Also, Larceny is rumoured to be very fast. Both can compile dynamically generated code into machine code on the fly, which is useful for GP. [TODO: show how to set up a real GP experiment with Ikarus]
The next step is to allow migration between colonies running on different computers or CPUs. Work in progress - coming soon.
I'm working on some macros to make fitness functions easy to construct.
You can't just use eval directly, since (1) it's really slow,
which is a problem for typical GP problems where you may want to test your
individuals in some kind of iterative process and
(2) it doesn't give you the right kind of control over the environment - you
typically want to provide functions such as %.
; macro allowing for lambda construction in custom environment
(define-syntax eval-with-bindings
(syntax-rules ()
((_ ((name value) ...) expression)
((eval `(lambda (name ...) ,expression)
(environment '(rnrs)))
value ...))))
(define (my-fitness-function individual)
; First we turn our s-expression into a procedure - this is where
; we decide what to provide in its environment, and also which parameters
; it takes. We can also provide macros to the function, by using let-syntax
; around the lambda expression.
(define function (eval-with-bindings ([% my-safe-divide]])
[foo 1024])
`(lambda (x y)
,individual))
; Now we can do whatever we need to do to evaluate 'function', which will
; be screaming fast native machine code under Ikarus or Larceny.
...)
This allows for the individual to be passed in as an s-expression, which is converted into a function which is compiled to machine code (if you're using Ikarus) in your hand-constructed environment. (Thanks to the guys on the ikarus-users list for helping me understand how to do this).
How can I get rid of the mutation (set!) in the functions
tree-ref and tree-combine?