Genetic Programming

<< "GPinit`"

Robert B. Nachbar

Evolutionary computing has become a popular method to robustly  search very large solution spaces for optimal or near-optimal solutions to a wide variety of problems. These methods begin with an initial, usually randomly generated, population of individual instances of a data structure that can be evaluated to a solution for the problem at hand. The term ‘solution’ is used rather loosely here because these primordial individuals are generally not very good solutions. Then, in the spirit of survival of the fittest, they are propagated directly, through recombination with mates, or by mutation into new, and hopefully more fit, individuals. This iterative process is repeated until one either grows weary or a satisfactory solution is found. The governing principle of these stochastic methods is that the fittest individuals possess a part of the final solution and by operating on them to make more-fit individuals one will eventually discover the optimal solution. Unlike the more traditional analytic methods (e.g., least squares, conjugate gradients, etc.), the stochastic nature of the multiple search trajectories endows these methods with the ability to avoid local minima and to sample a large portion of the solution space.

John R. Koza has recently described his approach to genetic programming with the Lisp language [Koza 1992]. It is difficult to provide a more thorough introduction to this topic than his, and the reader is directed to this work for further examples and more detailed discussions. What the present author intends to show here is that genetic programming is easily done in Mathematica as well.

Evolutionary Computing

In a genetic algorithm (GA) the individuals are character strings in whose elements are encoded the parameters that specify a solution for the problem at hand. A very close analogy between the strings and chromosomes has been drawn [Goldberg 1989]. In a GA, mutation is achieved by randomly selecting one of the letters of a string and exchanging it for a different letter of the alphabet. Typically a binary alphabet is used (i.e., 0 and 1). Recombination via mating is carried out by crossing over part of one parent with the homologous part of the other parent. This is done by selecting a point at random along the string and then exchanging the distal portions between the parents, thus generating two new offspring. Because the string is a coded representation of a solution, its length is fixed. In spite of this rigid linear data structure, GAs have been applied successfully to a large number of problems. See [Goldberg 1989] and [Davis 1991] for further details and examples. James Freeman has recently described the use of Mathematica for GAs [Freeman 1993].

A genetic program (GP) shares most of the features of a GA. However, the data structure that GP uses is hierarchical rather than linear. In addition, instead of employing an alphabet to encode the solution parameters, the parameters themselves are stored in the data structure along with the functions that operate on them. A graph-theoretical tree is suitable hierarchical data structure for a GP. (This is the same data structure that most compilers use to parse mathematical expressions.) The parameters occupy the terminal or external nodes (sometimes called leaves) and the functions occupy the internal nodes (of which, the highest is known as the root). These trees are customerily drawn "upside-down" with the root at the top and the leaves at the bottom. For example, the expression 2 sin x + b x + c can be represented as the tree shown in Figure 1.

DrawExpression[plus[times[2, sin[x]], times[b, x], c]];

This expression tree is in fact a reusable program whose inputs are the symbolic terminal nodes and whose output is the result produced by the root node.

In an expression tree, a subtree is the analog of a substring in a genetic algorithm. The genetic operations of mutation and crossover are applied to the subtrees.

Mathematica is a functional language that employs a hierarchical data structure for mathematical expressions. In addition, the language allows one to freely manipulate parts of expressions. These features make Mathematica eminently suited for genetic programming.

A Problem

We first became interested in genetic programming for solving symbolic regression problems, a generalization of curve-fitting. In traditional methods, such as least squares, one uses a presumed functional form (e.g., y  = mx + b) and determines only the coefficients (m and b) that best fit the observed data. In symbolic regression, one is interested in finding the best functional form as well as the coefficients. For example, given the series 5, 31, 121, 341, 781, 1555, 2801, 4681, what is the next entry? What is the 12th? To answer these questions, one must discover the the relationship between i, the position in the series, and Si, the i-th entry. The task of genetic programming is to find the appropriate combination of functions, variables, and coefficients to exactly reproduce this series. Then we can use the result to find any element of the series.

Expression Trees

Functions and Terminals

The basic data structure of a genetic program is the expression tree. At the root of the tree is a function, and the pendant nodes are the operands or arguments of the function. The internal representation of expressions in Mathematica (as revealed with the FullForm function) is just this sort of data structure.

FullForm[(x + 1) (y - z)]

A powerful feature of the tree structure is that it is recursive. Each subtree is itself a valid tree. All the internal nodes (Times, Plus, Plus, and Times in this instance) are at the roots of subtrees. Even the terminal nodes (1, x, y, -1, and z) are themselves trivial expression trees. From this example, it is obvious that functions occupy internal nodes and numbers and variables occupy terminal nodes.

Evaluation of an expression tree is carried out recursively. For a good description of trees and recursive programming see [Gaylord et al. 1993]. Since we are dealing with Mathematica expressions directly, we will not need a special evaluation function; Mathematica will do it automatically.

Notice that in the example above, the expression has been rearranged by Mathematica into a canonical form equivalent to (1 + x) * (y + (-1 * z)). The input form that matches the original expression is Times[Plus[x, 1], Subtract[y, z]]. This rearrangement is relatively minor and could be ignored. However, Times and Divide present difficulties because their use can introduce a new function, namely Power.

FullForm[x x]

FullForm[(x/y)]

This may be undesirable for the problem at hand because genetic recombination of individuals could easily replace the second argument of Power with an expression that does not evaluate to an integer. Aside from this, it introduces a general lack of control over the function set.

Another difficulty arises when we use built-in functions whose arguments are entirely made up of numbers. We loose the full tree structure and are left with but a single node.

FullForm[5 (3 + 2)]

One way of overcoming these problems is to keep the functions and their arguments in lists, as described by [Gaylord et al. 1993] (e.g., {Times, {5}, {Plus, {3}, {2}}}). Alternatively, we can introduce our own arithmetic functions (e.g., plus, times, subtract, and divide) that do not evalute to anything. Both methods preserve the full structure of the expression and both require a function that can evaluate the expression and return a usable result. The evaluator for the former representation must recursively traverse the nested lists while it builds up a Mathematica expression. For the latter method the evaluator need only provide a set of replacement rules. We have found the second approach more convenient, and it probably requires less memory as well.

evalRules = {plus -> Plus, times -> Times, 
        subtract -> Subtract, divide -> Divide}; Eval[expr_] := expr /. evalRules

Thus, for the examples above:

times[plus[x, 1], subtract[y, z]]

Eval[%]

times[x, x]

Eval[%]

divide[x, y]

Eval[%]

times[5, plus[3, 2]]

Eval[%]

At this point it is appropriate to discuss two requirements of the sets of functions and terminals used in GP: closure and sufficiency. The value at any terminal node or result of any function at an internal node can be an argument to the function at the next level up the tree. One must therefore ensure that the choice of terminals and functions be fully compatible so that functions can accept any arguments they may receive and still yield valid results. That is, the functions should be well defined and closed. For example, one would not want to freely mix numerical and Boolean terminals, as the result is still partially unevaluated:

FullForm[Eval[plus[2, True]]]

Nor does one want to encounter undefined results:

FullForm[Eval[divide[subtract[x, x], 0]]]

This latter situation (both the division by 0 and the resulting 0*Infinity) is easily remedied by using protected division which tests the denominator before the division is carried out. If it is zero, then a value such as 1 can be returned for the result which maintains closure [Koza 1992, p. 82]. However, it is more practical to return a very large machine number because it not only satisfies closure but is more appropriate in real-world problems than is 1 as it is more like the true result (∞) [Lee 1994]. (We could allow 1/0 to evaluate to Infinity and Infinity/Infinity to evalute to Infinity, but this causes problems further down the road if we try to compile the result.) Similar protected functions should be used for Sqrt and Log, and Exp should be protected from over- and underflow.

bigInteger = (2^63) - 1; bigReal = N[bigInteger]; PDivide[n_Integer, 0] := bigIntegerPDivide[n_Integer, 0. ] := bigIntegerPDivide[n_, 0] := bigIntegerPDivide[n_Real, 0] := bigRealPDivide[n_Real, 0. ] := bigRealPDivide[n_, 0. ] := bigRealPDivide[n_, d_?NumberQ] := (n/d)

evalRules = 
      ReplacePart[evalRules, PDivide, Position[evalRules, Divide]]

Eval[divide[subtract[x, x], 0]]

All of the type-testing for PDivide is necessary so that actual division does not take place before it is known if the denominator is 0 and that an appropriate precision result is returned.

Eval[divide[subtract[x, 1], y]]

Sufficiency of the function set and terminal set means that some combination of them is capable of producing an expression that is the solution to the problem. This responsibility lies solely with the user and is not always a straight forward task. If variables without sufficient explanatory power are used, the solution cannot be found. Likewise for functions. One must also be careful not to include extraneous functions or terminals as this can hamper the performance of the search for the solution.

Conditional Functions

Logical tests are very common in procedural programs written in languages such as FORTRAN or C, but somewhat unexpected in mathematical expressions. The If function in Mathematica not only controls which branch is followed, but it also returns the final value of the branch followed. Because the first argument of If must be a Boolean we cannot use it directly, but the functionality we seek can be implemented by using a modified If function that uses only numerical arguments. For example:

Attributes[ltz] = {HoldRest}; AppendTo[evalRules, 
      ltz[test_, t_, f_] :> If[test < 0, t, f]];

Eval[ltz[-3, x, y]]

Eval[ltz[2, x, y]]

The attribute HoldRest is necessary to prevent Mathematica from evaluating t and f prior to comparing test against 0. The reason is not obvious from this example, but if either t or f produced a side effect that was used elsewhere, one would want only the correct side effect to be produced. (The reader might try clearing the attribute and including a Print in each of the branches to see the effect.)

Constant Terminals

Frequently one needs constants to fully describe the solution to a problem. Even if one does not include them in the terminal set, they can be constructed during the evolution of the run. Integer and rational constants are the easiest to come across.

Eval[plus[x, x]]

Irrational constants can also be spontaneously generated. In a GP used to find trigonometric identities [Koza 1992, pp. 242-5], the constant π/2 was approximately constructed as

N[Eval[2 - Sin[Sin[Sin[Sin[Sin[Sin[Sin[Sin[1]] Sin[Sin[1]]]]]]]]]]

N[(π/2)]

To facilitate the incorporation of constants, the ephemeral random constant was introduced. It is used as a terminal during the creation of the initial generation, and whenever it is encountered, a random number of the appropriate type is generated. During the evolution of the population, these constants can be recombined in many ways to form new constants.

ephemeralReal := Random[Real, {-1, 1}]

Table[ephemeralReal, {5}]

Algebraic Simplification

As mentioned above, we chose to define our own arithmetic functions so that we could maintain better control of the functions used. With this change we also lost all of Mathematica’s built in simplification. For example,

(x - x) + 2 ((x + 2)/2)

whereas,

plus[subtract[x, x], times[2, divide[plus[x, 2], 2]]]

Just how much simplification is necessary is not always clear. The result 2 + x is certainly concise, which is desirable for a solution to a regression problem. However, we have also lost a great deal of genetic flexibility. For example, there are only 3 positions available at which genetic operations may occur in the simplified result as opposed to 11 in the original expression.

We will hope to strike a useful balance by just employing the associative,  identity, and inverse properties of algebra. Associativity is conferred by flattening out immediate subexpressions with the same head. We could have given the attribute Flat to the functions, but then plus[x] would not reduce to x.

plus[times[u, x], plus[y, z]]

plus[a___, b_plus, c___] := 
      Flatten[Unevaluated[plus[a, b, c]], 1, plus]plus[a_] := atimes[a___, b_times, c___] := 
      Flatten[Unevaluated[times[a, b, c]], 1, times]times[a_] := a

plus[times[u, x], plus[y, z]]

Very often we have seen expressions that contain subexpressions that will always reduce to 0 or 1 (e.g., subtract[x, x]). It is probably useful to provide rules for these cases. The following evaluation rules below make use of the identity and inverse properties of addition and multiplication:

plus[a___, 0, b___] := plus[a, b]plus[a___, 0. , b___] := plus[a, b]subtract[a_, 0] := asubtract[a_, 0. ] := asubtract[a_, a_] := 0times[a___, 1, b___] := times[a, b]times[a___, 1. , b___] := times[a, b]divide[a_, 1] := adivide[a_, 1. ] := adivide[a_, a_] := 1

Finally, here is the full effect:

plus[subtract[x, x], times[2, divide[plus[x, 2], 2]]]

Adam & Eve

The initial population in a GA is created by generating random character strings of a prescribed length, which is a fairly trivial task. In GP, on the other hand, we must generate random expression trees which have not only a length (depth), but also breadth. The inputs to a random expression generator are the lists of functions and terminals to use. In addition, the number of arguments that each function takes is needed. (For simplicity, we restrict plus and times to two arguments even though Mathematica has no such restriction for Plus and Times.)

funcs = {{plus, 2}, {subtract, 2}, {times, 2}, {divide, 2}}; terms := {x, y, z, Random[Integer, 3]};

We use recursion to construct the expression tree. The depth of the tree is controlled by decrementing a counter as we enter each level. The attribute HoldRest is necessary so that Random can be used in the terminal list.

randomElement[list_] := 
      list[[Random[Integer, {1, Length[list]}]
          ]]

RandomExpression[max_, funcs_, terms_] := 
      Module[{f, n}, {f, n} = randomElement[funcs]; 
        If[n > 0, f@@Table[RandomExpression[max - 1, funcs, terms], {n}], f]]
        RandomExpression[0, funcs_, terms_] := randomElement[terms]Attributes[RandomExpression] = {HoldRest};

RandomExpression[2, funcs, terms]

All the expressions produced in this manner have the same shape. They are full trees, that is, the distance from the root to each leaf is the same, and they always have a maximal number of nodes (that is, of course, before simplification). This was accomplished by restricting the use of terminals to the last level. Variably shaped trees are produced when terminals are allowed at any level.

comb = Join[funcs, ({#1, 0}&)/@terms];

re = RandomExpression[3, comb, terms]

re = times[plus[x, y], divide[times[y, x], plus[1, x]]];

The following two recursive functions allow us to ascertain the variability of the trees in the constructed population. This is important because the early use of terminals may prevent even one branch of a tree from reaching the prescribed depth.

DepthExpression[_[args__]] := 1 + Max[DepthExpression/@{args}]DepthExpression[_?AtomQ] := 0

SizeExpression[_[args__]] := 1 + Plus@@SizeExpression/@{args}SizeExpression[_?AtomQ] := 1

SizeExpression[re]

In constructing the initial population of expression trees, one should strive for variety. In the ramped-half-and-half method [Koza 1992, p. 93], the population is divided into equal groups for each depth and half the expression trees in each group are full and the other half are not. The function makePop does this. The attribute HoldAll is necessary so that Random can be used in the terminal set.

makePop[funcs_, terms_, nPop_, maxDepth_, minDepth_:  2] := 
      Module[{comb = Join[funcs, ({#1, 0}&)/@terms], pop, 
          d = maxDepth - minDepth + 1, i, depth, r}, 
        Table[depth = minDepth + Floor[(((i - 1) d)/nPop)]; 
          If[OddQ[i], 
            While[DepthExpression[
                  r = RandomExpression[depth, funcs, terms]] < depth]; r, 
            While[DepthExpression[r = RandomExpression[depth, comb, terms]] < 
                depth]; r], {i, nPop}]]Attributes[makePop] = {HoldAll};

pop = makePop[funcs, terms, 500, 6];

We can use DepthExpression and SizeExpression to check our generative method to make sure that we achieved the variety we sought. For a depth of 6 and all the functions taking 2 arguments, there will be a maximum size of 127 nodes (i.e., 27-1). First we load the standard package Statistics`DataManipulation` so that we can use the functions CategoryCounts and RangeCounts.

Needs["Statistics`DataManipulation`"];

depths = DepthExpression/@pop; d = Range[Min[depths], Max[depths]]; c = CategoryCounts[depths, d]; TableForm[Transpose[{d, c}], 
      TableHeadings -> {None, {"depth", "# trees"}}, 
      TableSpacing -> {0, 3}]

sizes = SizeExpression/@pop; s = Range[0, 128, 8]; c = RangeCounts[sizes, s]; s = Transpose[{Drop[s, -1], Drop[s, 1]}]; s = Apply[ToString[#1]<>"<=size<"<>ToString[#2]&, s, {1}]; TableForm[Transpose[{s, Take[c, {2, Length[c] - 1}]}], 
      TableHeadings -> {None, {"size", "# trees"}}, 
      TableSpacing -> {0, 3}]

Cain & Abel

The evolution of the population is carried out by applying the genetic operators to selected individuals to create offspring that become the next generation. In terms of frequency of occurance, reproduction and crossover are the two main genetic operators in GP, with mutation occuring less often. Reproduction, as the name suggests, is a direct copying of an individual expression from one generation to the next, crossover exchanges subexpressions from two parents to create two children, and mutation replaces a subexpression in an individual with a randomly generated subexpression. We also use a fourth operator, constant perturbation [Spencer, 1994], in those applications that employ constants as part of the terminal set.

Crossover

Crossover in GP is a bit more complicated than in a GA. To carry it out, two parents are selected at random from the population. Then a subexpression is randomly selected from one parent and exchanged with a randomly selected subexpression from the other parent.

Given two expression trees drawn on paper, a pair of scissors, and some glue, a person has very little trouble performing a crossover. However, this was the trickiest function to program for our GP. The process can be divided into two parts: first identifying the subexpressions to exchange, and second swapping them. We can use the built-in function Position to return a list of positions of subexpressions which are the parts that match the very simple pattern _.

Adam = RandomExpression[2, funcs, terms]

Adam = plus[divide[z, x], times[y, 2]];

p = Position[Adam, _, Heads -> False]

Next we use the built-in function Part to excise a subexpression. It takes as its second and following arguments the indices of the part, so all we have to do is make a list of the indices of all the subexpressions in the expression and then select one of them at random. To use one of the lists of indices returned by Position, it has to be recast as a sequence.

Adam
      [[Sequence@@
        p[[6]]]]

The built-in  function ReplacePart can be used to insert the excised subexpressions. It takes as its third argument the position for replacement (we have to make a slight amendment to it so that the whole expression can be replaced). Finally, we present below the full crossover function.

Unprotect[ReplacePart]; ReplacePart[expr_, new_, {}] := newProtect[ReplacePart];

CrossOver[parents_] := 
      Module[{pos, ind, sub}, 
        pos = (Position[#1, _, Heads -> False]&)/@parents; 
        ind = randomElement/@pos; 
        sub = MapThread[#1
                [[Sequence@@#2]]&, {
              parents, ind}]; 
        MapThread[ReplacePart, {parents, Reverse[sub], ind}]]

Now to show it in action, let's create Adam's mate, Eve, and then produce their offspring, Cain and Abel.

Eve = RandomExpression[3, funcs, terms]

Eve = 
      divide[subtract[plus[1, 1], times[y, x]], 
        divide[subtract[z, x], plus[z, 1]]];

{Cain, Abel} = CrossOver[{Adam, Eve}];

Cain

Cain = plus[divide[z, x], x];

Abel

Abel = 
      divide[subtract[plus[1, 1], times[y, x]], 
        divide[subtract[z, times[y, 2]], plus[z, 1]]];

As one can see in Figure 2, part {2} of Adam (times[y, 2]) was swapped with part {2, 1, 2} of Eve (x).

Needs["Graphics`Arrow`"]; Show[
      GraphicsArray[{{
            DrawExpression[Adam, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{2, Red}}], 
            DrawExpression[Eve, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{2, 1, 2}, Blue}}]}, {{}, {}}, {
            DrawExpression[Cain, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{2, Blue}}], 
            DrawExpression[Abel, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{2, 1, 2}, Red}}]}}, 
        Epilog -> {
            Line[{{0.5, 0.90}, {0.5, 0.83}, {1.5, 0.83}, {1.5, 0.90}}], 
            Line[{{1.0, 0.83}, {1.0, 0.62}}], 
            Line[{{0.5, 0.62}, {1.5, 0.62}}], 
            Arrow[{0.5, 0.62}, {0.5, 0.5}, HeadScaling -> Absolute], 
            Arrow[{1.5, 0.62}, {1.5, 0.5}, HeadScaling -> Absolute], 
            Text["Adam", {0.45, 0.90}, {1, 0}], 
            Text["Eve", {1.55, 0.90}, {-1, 0}], 
            Text["Cain", {0.45, 0.5}, {1, 0}], 
            Text["Abel", {1.55, 0.5}, {-1, 0}]}]];

There are some interesting special cases that should be discussed. If the exchange takes place between two terminal nodes, then we have effected two mutations. Even if crossover takes place between two identical parents, the offspring will in general not be identical because the subexpressions selected for exchange will not necessarily be the same. This is in contrast to a GA where the fixed data structure forces the generation of two offspring not only identical to each other, but also identical to the parents. During the evolution of the population, the depth of individuals tends to grow because subexpressions of different depths are exchanged. While this will allow great flexibility in finding a solution, these larger expressions take longer to evaluate. Therefore it is practical to impose an upper limit on the depth of an expression. In the event that a crossover will exceed this limit, then the offending offspring is replaced with one of its parents [Koza 1992, p. 104], as shown in this revised definition of CrossOver:

Clear[CrossOver]; CrossOver[parents_, maxDepth_:  16] := 
      Module[{pos, ind, sub, children}, 
        pos = (Position[#1, _, Heads -> False]&)/@parents; 
        ind = randomElement/@pos; 
        sub = MapThread[#1
                [[Sequence@@#2]]&, {
              parents, ind}]; 
        children = MapThread[ReplacePart, {parents, Reverse[sub], ind}]; 
        MapThread[
          If[DepthExpression[#1] ≤ maxDepth, #1, #2]&, {children, 
            parents}]]

Mutation

In the traditional GA, mutation merely exhanges the letter in a randomly selected gene for another one from the allowed alphabet. Typically, a bit is flipped, as most GAs use a binary representation. In genetic programming, mutation replaces a randomly selected subexpression with a new, randomly generated subexpression. The depth of the new subexpression can be anywhere from 0 (i.e., a terminal) to some maximum based on the overall maximum depth for an expression, which may be larger than the orignal subexpression.

Mutate[expr_, funcs_, terms_, maxDepth_:  16] := 
      Module[{pos, ind, depth, newSub}, 
        pos = Position[expr, _, Heads -> False]; 
        ind = randomElement[pos]; 
        depth = Random[Integer, (1/4) (maxDepth - Length[ind])]; 
        newSub = RandomExpression[depth, funcs, terms]; 
        ReplacePart[expr, newSub, ind]]Attributes[Mutate] = {HoldRest};

In the example shown in Figure 3, we see that part {2,1} of Simon of depth 0 was replaced by a subexpression of depth 1 to become Peter.

Simon = RandomExpression[2, funcs, terms]

Simon = times[divide[y, z], subtract[1, z]];

Peter = Mutate[parent, funcs, terms]

Peter = times[divide[y, z], subtract[times[z, y], z]];

Show[
      GraphicsArray[{{
            DrawExpression[Simon, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{2, 1}, Red}}]}, {}, {
            DrawExpression[Peter, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{2, 1}, Red}}]}}, 
        Epilog -> {
            Arrow[{0.5, 0.825}, {0.5, 0.5}, HeadScaling -> Absolute], 
            Text["Simon", {0.45, 0.825}, {1, 0}], 
            Text["Peter", {0.55, 0.5}, {-1, 0}]}]];

Constant Perturbation

In many instances, an expression tree can be improved by making a small change to one of the component constants. Fine tuning of this sort is often more efficient than relying on mutation to achieve the same effect. We begin by identifying only those terminal parts of an expression tree that hold a constant numeric value with the pattern _?NumberQ.

terms := {x, Random[Integer, {-5, 5}]};

Saul = RandomExpression[3, funcs, terms]

Saul = 
      divide[times[subtract[4, x], times[3, -1]], 
        subtract[subtract[-5, x], subtract[3, x]]];

cp = Position[Saul, _?NumberQ]

(Saul[[Sequence@@#1]]&)/@cp

Now one of these nodes is randomly selected and the value stored there is perturbed by a small amount. Note that we replace Reals with Reals and Integers with Integers. Similar actions should be taken for Complex and Rational constants as well.

perturbConstant[n_Real] := Random[Real, n {0.8, 1.25}]perturbConstant[n_Integer] := 
      Random[Integer, {Min[n - 3, Floor[0.8 n]], 
            Max[Ceiling[1.25 n], n + 3]}] /; n ≥ 0perturbConstant[n_Integer] := 
      Random[Integer, {Min[n - 3, Floor[1.25 n]], 
            Max[Ceiling[0.8 n], n + 3]}] /; n < 0

Perturb[expr_] := 
      Module[{constants, ind, oldTerm, newTerm}, 
        constants = Position[expr, _?NumberQ]; 
        If[Length[constants] > 0, ind = randomElement[constants]; 
          oldTerm = 
            expr[[Sequence@@ind]]; 
          newTerm = perturbConstant[oldTerm]; 
          ReplacePart[expr, newTerm, ind], expr]]

The example in Figure 4 shows that the value 4 at part {1, 1, 1} of Saul was replaced with the value 1 in Paul.

Paul = Perturb[Saul]

Paul = 
      divide[times[subtract[1, x], 3, -1], 
        subtract[subtract[-5, x], subtract[3, x]]];

Show[
      GraphicsArray[{{
            DrawExpression[Saul, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{1, 1, 1}, Red}}]}, {}, {
            DrawExpression[Paul, DisplayFunction -> Identity, 
              xMargin -> 1.2, deltaX -> 1.5, 
              style -> {{{1, 1, 1}, Red}}]}}, 
        Epilog -> {
            Arrow[{0.5, 0.45}, {0.5, 0.25}, HeadScaling -> Absolute], 
            Text["Saul", {0.47, 0.45}, {1, 0}], 
            Text["Paul", {0.53, 0.25}, {-1, 0}]}]];

Fitness & Selection

How does one judge the usefulness of an expression tree for solving a problem? In the problem outlined above, we have a set of observations or cases against which we can measure the performance of the program.

cases = {{1, 5}, {2, 31}, {3, 121}, {4, 341}, {5, 781}, {6, 1555}, {
          7, 2801}, {8, 4681}}; x = (Drop[#1, -1]&)/@cases; y = Last/@cases;

The better the program does, the more fit it is and thus likely to survive to the next generation and/or mate with other individuals. What constitutes ‘better’ depends on the problem, but most commonly, we want to minimize an error function. In this case the error is the raw fitness. One could use the sum of absolute errors or the sum of squared errors over the fitness cases, and there is no compelling reason to choose one over the other. In the functions SumAbs and SumSqr, the arguments x and y are the independent and dependent values of the fitness cases, vars is the list of variables or terminals that correspond to the columns in x, and expr is an expression.

SumAbs[x_, y_, vars_, expr_] := 
      Module[{f}, f = Function[vars, Evaluate[Eval[expr]]]; 
        N[Plus@@Abs[y - Apply[f, x, {1}]]]]

SumSqr[x_, y_, vars_, expr_] := 
      Module[{f, d}, f = Function[vars, Evaluate[Eval[expr]]]; 
        d = N[y - Apply[f, x, {1}]]; d . d]

funcs = {{plus, 2}, {subtract, 2}, {times, 2}, {divide, 2}}; terms = {i, 1, 2, 3, 4}; vars = {i};

pop = makePop[funcs, terms, 10, 6];

pop = {times[subtract[3, 4], plus[i, 3]], 
        times[divide[i, 2], divide[4, 3]], 
        plus[times[plus[3, 3], divide[4, 2]], 
          plus[divide[i, i], times[2, 4]]], 
        subtract[divide[3, plus[4, i]], divide[divide[i, 4], 2]], 
        plus[subtract[subtract[plus[3, 1], times[i, 2]], 
            divide[plus[2, 4], divide[3, 3]]], 
          subtract[divide[divide[4, i], plus[2, 3]], 
            plus[times[4, 1], times[1, i]]]], 
        times[divide[subtract[4, 2], 3], times[times[divide[1, 4], 1], 1]], 
        divide[subtract[
            times[subtract[times[i, 4], divide[4, i]], 
              times[plus[4, 4], subtract[1, 3]]], 
            divide[times[subtract[4, 3], divide[2, 2]], 
              subtract[divide[i, 2], subtract[i, 2]]]], 
          subtract[
            times[times[divide[1, i], subtract[3, 4]], 
              times[subtract[4, 2], subtract[i, 4]]], 
            subtract[divide[divide[3, 3], subtract[i, 1]], 
              plus[divide[3, 1], subtract[i, i]]]]], 
        plus[times[times[i, 2], 
            plus[subtract[2, plus[1, 2]], times[3, plus[3, i]]]], 3], 
        plus[divide[
            subtract[
              subtract[subtract[subtract[1, 2], times[1, i]], 
                subtract[divide[i, i], plus[2, i]]], 
              subtract[plus[subtract[i, 2], subtract[i, 1]], 
                divide[times[i, 3], subtract[2, 3]]]], 
            divide[times[times[times[4, 1], subtract[2, i]], 
                subtract[divide[1, i], times[i, i]]], 
              times[times[plus[1, 2], plus[2, 4]], 
                subtract[subtract[4, 2], plus[3, 4]]]]], 
          subtract[
            plus[times[times[divide[2, 3], divide[4, i]], 
                times[subtract[3, i], divide[1, i]]], 
              plus[divide[plus[2, 2], divide[2, 4]], 
                times[times[4, 3], subtract[1, 3]]]], 
            divide[divide[divide[divide[1, i], plus[1, 4]], 
                plus[plus[4, 1], divide[4, 2]]], 
              plus[plus[plus[2, 4], subtract[4, 3]], 
                times[subtract[i, 4], divide[i, i]]]]]], 
        times[4, 
          subtract[i, 
            times[plus[subtract[4, i], 2], 
              subtract[1, divide[times[1, 2], 2]]]]]};

The list of raw fitnesses is obtained by mapping sumAbs onto the population.

rawFit = (SumAbs[x, y, vars, #1]&)/@pop

One should be aware that it is always possible to minimize error by adding more adjustable parameters to the fitting function, and we have done nothing to protect against this so-called overfitting. One could build parsimony into the fitness function by inflating the error proportional to the number of terminals and/or functions in the expression tree, however, it has been shown [Koza 1992, pp. 612-4] that penalizing the fitness by the expression size seriously degrades performance. Therefore, it is more convenient to enforce parsimony via the genetic operators and algebraic simplification, as described above.

Depending on the problem, we may want to either minimize or maximize the raw fitness. To simplify things, the raw fitness is converted to standardized fitness, which is always minimized. When raw fitness is error, then this definition suffices.

standardizedFit[rawfit_] := rawfit

standFit = standardizedFit[rawFit]

When the raw fitness is a score, as in a game, which should be maximized, then the following definitions may be appropriate. If there is no theoretical upper bound to the score, maxScore may be set to some practical value or simply the maximum value of the population's raw fitness.

standardizedFit[rawfit_] := maxScore - rawfitstandardizedFit[rawfit_] := Max[rawfit] - rawfit

Finally, the standardized fitnesses are inverted to give the adjusted fitness.

adjustedFit[standfit_] := (1/(1 + standfit))

adjFit = adjustedFit[standFit]

The adjusted fitness is used to determine the likelihood an individual is selected for input to a genetic operator. The following example shows that it is not necessary to convert adjusted fitness to normalized fitness (Koza 1992, pp. 97-8). Adjusted fitness is used directly by the function selectOne (defined below). The last column in the table shows that the probability of selection is indeed the same as the normalized fitness.

af = Table[Random[], {naf = 20}]; cumSum = Rest[FoldList[Plus, 0, af]]; t = Table[selectOne[cumSum], {n = 10000}]; c = N[CategoryCounts[t, Range[naf]]/n]; TableForm[Transpose[{af, (af/Last)[cumSum], c}], 
      TableHeadings 
        -> {
          None, {"adjusted\nfitness", "normalized\nfitness", 
            "probability\nselected"}}, TableSpacing -> {0, 3}]

There are three main types of selection. In fitness proportionate selection the probability of selection is directly proportional to the relative fitness. This method is also known as roulette wheel selection. It is most easily accomplished by using a list of cumulative sums of adjusted fitnesses [Freeman 1993].

cumSum = Rest[FoldList[Plus, 0, adjFit]]

Selection is then done by generating a random number between 0 and the final cumulative sum, and finding the position of first cumulant that exceeds it. selectOne returns the index of the selected individual.

selectOne[list_] := 
      Module[{r = Random[] Last[list]}, 
        Position[list, _?(#1 ≥ r&), {1}, 1]
            [[1, 1]]]

selectOne[cumSum]

pop[[%]]

While this definition of selectOne is an example of functional programming, it is not very fast because each element of the list of cumulative fitnesses must be examined in turn until one matches the pattern. The algorithm is of order O(n). Because the list is nondecreasing, we can use a more efficient binary search algorithm of order O(log2 n). The procedural binary search version shown below is more than 3 times faster for a list of 500 elements.

fit = Table[Random[], {500}]; cumSum = Rest[FoldList[Plus, 0, fit]];

rs = $RandomState; t1 = First[Timing[Table[selectOne[cumSum], {100}]]]

Clear[selectOne]; selectOne[list_] := 
      Module[{lo = 1, hi = Length[list], mid, r = Random[] Last[list]}, 
        While[lo ≠ hi - 1, mid = Round[((lo + hi)/2)]; 
          If[r < list[[mid]], hi = mid, 
            lo = mid]]; 
        If[r ≤ list[[lo]], 
          lo, hi]]

$RandomState = rs; t2 = First[Timing[Table[selectOne[cumSum], {100}]]]

(t1/t2)

Rank selection is closely related to fitness proportionate selection. In this method, the probability of selection is proportional to only the rank, or order, of the adjusted fitness and not the actual magnitude. It would be used when one wants to enhance the distinction between individuals that are nearly equally fit. The function rankedFit takes a list of adjusted fitnesses and returns their ranks (where 1 is low and ties are allowed).

toRank[adjFit_] := 
      Module[{f = adjFit ⋃ (SameTest -> SameQ), c, 
          r = Range[Length[adjFit]], g, b, nr}, 
        c = CategoryCounts[f, adjFit]; 
        g = ((b = Take[r, #1]; r = Drop[r, #1]; b)&)/@c; 
        nr = (Plus@@#(1/Length)[#1]&)/@g; 
        nr[[Flatten[(Position[f, #1]&)/@adjFit]
            ]]]

N[toRank[adjFit]]

The cumulative sum of ranks would then  be used in place of the cumulative sum of adjusted fitnesses in selectOne.

cumSum = Rest[FoldList[Plus, 0, %]]

The third method of selection is tournament selection. Here, two or more individuals are drawn at random from the population and the most fit one is selected.

Putting it all together

Before we present the synthesis of geneticProgram, we need to define a success predicate that returns true when an optimal or, depending on the problem, near-optimal solution is found. When the raw fitness is error, then a perfect solution would have zero error, otherwise we might be willing to accept an error below a predefined threshold, as shown in this example:

successQ[rawfit_] := 
      Position[rawfit, _?(#1 < 0.5&)] ≠ {}

We would like GeneticProgram to be as generic as possible so that it can be used in many problem domains without the need of alteration. The number of required arguments should be kept to a minimum, so we have made liberal use of options.

Options[GeneticProgram] = {PopulationSize -> 500, 
        MaxGenerations -> 50, 
        GeneticOpers 
          -> {{CrossOver, 2,  .85}, {Reproduce, 1,  .01}, {Mutate, 
              1,  .01}, {Perturb, 1,  .13}}, MaxDepthInitial -> 6, 
        MaxDepthCreated -> 16, Selection -> Roulette, 
        FitnessFunction -> SumAbs};

\(GeneticProgram[fitCases_, vars_, funcs_, terms_, standardFit_, 
        successQ_, options___] := 
      Module[{popSize, maxGen, genOpers, cumProb, maxInitial, maxCreated, 
          optlist, result, pop, g, selection, x, y, rawFit, adjFit, cumSum, 
          p, newPop, i, j, k, np, op, parents, children, bGen, bRun}, 
        optlist = Flatten[{options}]; 
        popSize = \(PopulationSize /. optlist\) /. Options[GeneticProgram]; 
        maxGen = \(MaxGenerations /. optlist\) /. Options[GeneticProgram]; 
        maxInitial = 
          \(MaxDepthInitial /. optlist\) /. Options[GeneticProgram]; 
        maxCreated = 
          \(MaxDepthCreated /. optlist\) /. Options[GeneticProgram]; 
        genOpers = \(GeneticOpers /. optlist\) /. Options[GeneticProgram]; 
        cumProb = Rest[FoldList[Plus, 0, Apply[#3\/#2&, genOpers, {1}]]]; 
        genOpers = \((Drop[#1, \(-1\)]&)\)/@genOpers; 
        genOpers = 
          genOpers /. {Reproduce -> \((#1&)\), 
              Mutate -> 
                \(({Mutate[Sequence@@#1, funcs, terms, maxCreated]}&)\), 
              CrossOver -> \((CrossOver[#1, maxCreated]&)\), 
              Perturb -> \(({Perturb[Sequence@@#1]}&)\)}; 
        selection = \(Selection /. optlist\) /. Options[GeneticProgram]; 
        fit = \(FitnessFunction /. optlist\) /. Options[GeneticProgram]; 
        x = \((Drop[#1, \(-1\)]&)\)/@fitCases; y = Last/@fitCases; 
        pop = makePop[funcs, terms, popSize, maxInitial]; g = 0; 
        rawFit = \((fit[x, y, vars, #1]&)\)/@pop; 
        adjFit = \((1\/\(1 + #1\)&)\)/@standardFit[rawFit]; 
        p = \(Position[adjFit, Max[adjFit]]\)[[1, 1
            ]]; 
        bGen = {{rawFit[[p]], 
              pop[[p]]}}; 
        bRun = {adjFit[[p]], 
            rawFit[[p]], 
            pop[[p]]}; 
        Print["g = ", g, ", bestOfRun = ", 
          rawFit[[p]], 
          ", bestOfGen = ", 
          rawFit[[p]]]; 
        While[g < maxGen && \(! \((successQ[rawFit])\)\), 
          cumSum = 
            Which[selection === Roulette, Rest[FoldList[Plus, 0, adjFit]], 
              selection === Rank, Rest[FoldList[Plus, 0, toRank[adjFit]]]]; 
          newPop = Table[Null, {popSize}]; i = 0; 
          While[i < popSize, {op, np} = 
              genOpers[[selectOne[cumProb]
                  ]]; p = Table[selectOne[cumSum], {np}]; 
            parents = pop[[p]]; 
            children = op[parents]; 
            For[j = 1, j ≤ np && i < popSize, \(j++\), 
              newPop[[\(++i\)]] = 
                children[[j]]]; ]; 
          pop = newPop; \(g++\); rawFit = \((fit[x, y, vars, #1]&)\)/@pop; 
          adjFit = 1\/\(1 + standardFit[rawFit]\); 
          p = \(Position[adjFit, Max[adjFit]]\)[[1, 1
              ]]; 
          AppendTo[
            bGen, {rawFit[[p]], 
              pop[[p]]}]; 
          If[adjFit[[p]] > First[bRun], 
            bRun = {adjFit[[p]], 
                rawFit[[p]], 
                pop[[p]]}]; 
          If[Mod[g, 5] == 0, 
            Print["g = ", g, ", bestOfRun = ", 
              bRun[[2]], 
              ", bestOfGen = ", 
              rawFit[[p]]]]; ]; 
        If[Mod[g, 5] ≠ 0, 
          Print["g = ", g, ", bestOfRun = ", 
            bRun[[2]], 
            ", bestOfGen = ", 
            rawFit[[p]]]]; {
          BestOfRun -> Rest[bRun], BestOfGen -> bGen, 
          FinalFit -> rawFit, FinalPop -> pop}]\)Attributes[GeneticProgram] = {HoldAll};

Two interesting features of GeneticProgram should be pointed out. First, the genetic operators, their numbers of input expressions, and their relative frequencies are kept in a list, much the same way as the set of functions. This allows us to draw any one of them at random and apply it to the appropriate number of individuals drawn from the population. Second, the asexual genetic operators return a single expression, whereas the sexual one (crossover) returns a list of expressions. One of the great strengths of Mathematica is its ability to be customized on the fly. Thus we are able to recast Mutate and Perturb to return a list of one, and define Reproduce with the following snippet of code:

genOpers = 
      GenOpers /. {Reproduce -> (#1&), 
          Mutate -> 
            ({Mutate[Sequence@@#1, funcs, terms, maxCreated]}&), 
          CrossOver -> (CrossOver[#1, maxCreated]&), 
          Perturb -> ({Perturb[Sequence@@#1]}&)};

The relative frequencies for the genetic operators do not indicate the frequency with which they are invoked, but rather the frequency with which new individuals are created from the operations. Therefore, to obtain the necessary probability of using an operator, the relative frequency is divided by the number of parents (which is the same as the number of offspring produced).

The result returned by GeneticProgram is a list of replacement rules. BestOfRun is a couplet of the raw fitness and the expression for the best individual found during the entire course of the run. BestOfGen is a list of couplets of the raw fitness and the expression for the best individual found at each generation of the run. FinalFit is a list of the raw fitnesses for the final population, and FinalPop is the corresponding list of expressions. In many of the problems we encounter in medicinal chemistry, there is no one correct answer, so we prefer to look at an ensemble of good answers.

Trying It Out

Finding some solutions

So, back to the problem: finding an expression for our series 5, 31, 121, 341, 781, 1555, 2801, 4681. Here is the list of the eight observations:

cases = {{1, 5}, {2, 31}, {3, 121}, {4, 341}, {5, 781}, {6, 1555}, {
          7, 2801}, {8, 4681}};

The function set contains just the four arithmetic functions; nothing in the data warrants the use of trigonometric functions or logarithms.

funcs = {{plus, 2}, {subtract, 2}, {times, 2}, {divide, 2}}; terms := {i, Random[Integer, -3, 3]}

Standardized fitness is the same as raw fitness, and the success predicate will be a perfect fit.

standardizeFit[rawFit_] := rawFit

successQ[rawFit_] := Position[Chop[rawFit], 0] ≠ {}

We will use a population size of 500, a typical value. This number may seem rather high for such a simple problem, but it is required to allow the population to maintain the necessary diversity. It is difficult, however, to predict how the population size scales with the complexity of the problem. Suffice it to say here that we have solved a 3-variable problem with the same population size.

$RandomState = 
      2377645317229716664685496329582157415741425953655155128585616834948144146203580264738161970735009674346956811643336202116685336024086267847732568906324631266094917606622655511666474526420404120885261762711671400531973163288740165388202781377628915214942679702186428778647373571291602386191450906863331492720153415487265079100017556824688717951341153642937222582752608378742236121338627248338620004303284641604510435634460277032089206563259072802883540752277321687443969845309753573586295744218565378725918054793517262586732108526091819279299192809795573413387;

result = 
      GeneticProgram[cases, {i}, funcs, terms, standardizeFit, successQ, 
        FitnessFunction -> SumSqr, PopulationSize -> 500, 
        MaxGenerations -> 50];

Here is the best fitting expression and its depth:

b = BestOfRun /. result

DepthExpression[Last[b]]

DrawExpression[Last[b], AspectRatio -> (210/480)];

Here is the simplified form of the expression:

be = Eval[Last[b]]

bestFit = Together[Expand[be /. PDivide -> Divide]]

Now, let's look at the runners up. Sometimes the numerically best solution does not make good physical, chemical, or biological sense for the problem, so it is always a good idea to examine plausible alternatives. In the real-world where experimental error lurks, these second- and third-best answers may well be better than the best fit solution when one also takes into account parsimony, the physical interpretation of the GP-determined expressions, and test cases that were not used in training (i.e., validation).

fp = Sort[Transpose[{FinalFit, FinalPop} /. result]]; First/@Take[fp, 5]

Here is the second best fitting expression and its depth:

sb = fp[[2]]

DepthExpression[Last[sb]]

Here is the simplified form of the expression:

sbe = Eval[Last[sb]]

secondBestFit = Together[Expand[sbe /. PDivide -> Divide]]

Validating the solutions

Suppose now that we are given the next four entries in the series to test our results. (If one cannot obtain new data with which to validate ones result, then one should hold back a portion of the data from the fitting process for this purpose.) Thus, we have 7381, 11111, 16105, and 22621 for i = 9, 10, 11, and 12.

(bestFit /. i -> #1&)/@Range[9, 12]

Viola! Still a perfect fit. Here is a comparison of the two results:

TableForm[
      Transpose[{Range[12], (bestFit /. i -> #1&)/@Range[12], 
          (N[secondBestFit /. i -> #1]&)/@Range[12]}], 
      TableHeadings -> {None, {"i", "best", "second best"}}, 
      TableSpacing -> {0, 3}]

Closing Remarks

All the functions described here as well as several others that should be useful have been collected and made into the standard [Maeder 1991]  packages GeneticProgramming.m and SymbolicRegression.m. The options and output have been expanded to include information that we have found helpful in deciding whether or not to continue a run or to start over with a new population. We have attempted to isolate the problem specific functions and definitions in SymbolicRegression.m and to keep GeneticProgramming.m as generic as possible.

Currently, we exploring ways to evolve parsimonious expressions, employ smaller but more diverse populations, and to incorporate analytic methods to determine the values for constants.

The author wishes to acknowledge several helpful suggestions from the editor, and advice from Alan DeGuzman and Robby Villegas at the technical support desk at Wolfram Research, Inc.. He also wishes to thank Mary Curran Nachbar, a high school student, for careful reading of the manuscript to identify vague and imprecise terms and phrases.

References

Davis, L. 1991. Handbook of Genetic Algorithms. Van Nostrand Reinhold, New York.

Freeman, J. 1993. Simulating a Basic Genetic Algorithm. The Mathematica Journal3, 52-56.

Gaylord, R. J., Kamin, S. N., and Wellin, P. R. 1993. Introduction to Programming with Mathematica. Telos, Santa Clara, California.

Goldberg, D. E. 1989. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley, Reading, Massachusetts.

Koza, J. R. 1992.Genetic Programming: On the Programming of Computers by Means of Natural Selection. MIT Press, Cambridge, Massachusetts.

Lee, G. Y. 1994. private communication.

Maeder, R. 1991. Programming in Mathematica, Second Edition. Addison-Wesley, Reading, Massachusetts.

Spencer, G. 1994. Automatic Generation of Programs for Crawling and Walking. Advances in Genetic Programming,  Kenneth E. Kinnear, Jr., ed. MIT Press, Cambridge, Massachusetts.

Robert B. Nachbar
Merck Research Laboratories, R50S-100, P.O.Box 2000, Rahway, NJ 07065, (908)594-7795
nachbar@merck.com

The electronic supplement contains packages SymbolicRegression.m and GeneticProgramming.m and the notebook SymbolicRegression.ma. The latter uses the packages for the example in the previous section.