Genetic Programming with Mathematica

Background

(All of the work in this document is directly related to the bookGenetic Algorithms by John Koza. [Koza 1])

My interest in this topic began when I was reading Mastering Mathematica by John W. Gray, and I saw following quote (pg 169):

" 'Pascal is for building pyramids—imposing, breathtaking, static structures built by armies pushing heavy blocks into place. Lisp is for building organisms—imposing, breathtaking, dynamic structures built by squads fitting fluctuating myriads of simpler organisms into place.'
[Abelson].
						We, of course, intend to replace 'Pascal' by 'C' and 'Lisp' by 'Mathematica.' ".
 [Gray]
 
Before reading that, it had never occurred to me that Mathematica could be used for artificial intelligence. When I saw [Koza 1] in a bookstore, I knew its Lisp-based contents could somehow be converted to Mathematica.—I just didn’t realize it would come so easily. Keep in mind that I have absolutely no knowledge of Lisp, so all the following is based on the text of [Koza 1], but not the Lisp programming examples in the back.

Who this is for

This is written for anyone interested in implementing Genetic Algorithms with Mathematica.. You do not have to be an expert, but to understand this notebook you should read at least the first few chapters of [Koza 1]. If there are any other books that explain the purpose and theories behind Genetic Algorithms, they will also do. Most of the code in this document will apply to any problem involving Genetic Algorithms with little or no adaptation.

What is a Genetic Algorithm?

The purpose of this notebook is not to explain Genetic Programming, but if you want a brief explanation:
Genetic Algorithms are tools to solve problems. They provide the computer with a very intuitive, simple, and general way to solve problems from a very wide range of possibilities. There are not very many steps to know, and the major ones are:

1. Creating an initial random population.
2. Testing each member of the population for fitness.
3. Genetically breeding the most fit members to produce a new generation.
4. Repeating steps 2 and 3 until some sort of goal or limit is reached.

The intermediate steps are described in detail later, but these four are by far the most important to understand. This whole process is know as natural selection.

Implementation

Creating an initial population

When creating the initial population, diversity is the most important thing to encourage. It is relatively easy to prevent having duplicate members of the population, but that does not insure diversity (I actually did not include a checking function to ensure lack of duplicate members partly because I usually use smaller populations). It is very easy to have many individuals who are very similar, but different enough to pass an equality check with their peers. I have found that to have the most diversity, the individual–producing function must be able to create individuals with random depths and relationships of functions. There should be a consistent random element for creating every aspect of each individual.

Here is the basic function for creating a new population:

createRandomPop[maxLevels_,functionList_,varList_,numList_,popSize_,doPrint_]:=   {
	initPopulation = {},
	Do[ {If[doPrint == True, Print[i,":"] ],
		If[Random[Integer,{1,2}] == 1,
		levels = Random[Integer,{2,maxLevels}],
		levels = maxLevels],
		newIndividual = x,
	While[newIndividual == varList[[1]] || newIndividual == varList[[2]] || newIndividual == 1, 
		createRandomExp[levels,functionList,varList, numList, doPrint]
																	   ],
		initPopulation = Append[initPopulation, newIndividual]},
							{i, 1, popSize}]
												}

As you can see, this scheme uses global variables. I am aware of methods in Mathematica which can easily hide variables to other functions, but I figured that it was important to have functions that can easily exchange data and don't have prototypes a mile long. The variables I chose to include in the prototypes were those I thought would be necessary to use the function alone.
maxLevels is the maximum depth and individual can be. I decided to determine the depth in this function so it can be possible later to use the createRandomExp function to control the exact depth. The array functionList is a list of functions that can be used in the creation of each individual. The array numList is a list of the number of arguments allowed in each function. For example if your function list looked like:

functionList = {Plus, Abs}

then varList would be:

numList = {2, 1}

because addition normally requres two arguments and Abs only requires one. From the little I know of Lisp, I know that Mathematica functions internally are very similat to Lisp functions. In Lisp, Plus can only have two arguments. So 1 + 2 + 3 must be written as (+1 (+2 3)). By specifying "2" for the number of arguments for Plus, Mathematica is effectively simulating Lisp.

varList is an array of the variables (the actual inputs) that can be used in each individual. So if

varList = {x, y}

Then x and y are the two variables that are used for every individual created.

popSize is the number of individuals to be created.

doPrint is a boolean value that represents whether or not display each individual as it is created. Since each individual is stored in an array, a list can easily be displayed later.

Below are the actual functions used for creating new random individuals (referenced from createRandomPop). It is not important to know how they work, because the prototype of createRandomExp is very similar to that of createRandomPop. It is only important to know that they create a global named newIndividual (which is stored in an array by createRandomPop).

createRandomExp[levels_,functionList_,varList_,numList_,doPrint_]:=  {
	partNumber = Random[Integer,{1,Length[functionList]}],
	functionPart = functionList[[ partNumber ]],
	numberArgs = numList[[ partNumber ]], 
	argumentPart1 = varList[[ Random[Integer,{1,Length[varList]}] ]],
	argumentPart2 = varList[[ Random[Integer,{1,Length[varList]}] ]],
	If[numberArgs == 1,	fullFunction = functionPart[argumentPart1],
				fullFunction = functionPart[argumentPart1,argumentPart2]	
							],
	newIndividual = fullFunction,
		Do[ addPart[newIndividual,functionList,varList,numList]          
						,{i,1,(levels - 2)}]
					   If[doPrint == True, Print[newIndividual]]
							}

addPart[exp_,functionList_,varList_,numList_]:=  {
	getPos[exp];,
	thePoint = {};,
	Do[ { 
	thePoint = Append[thePoint,partList[i] ]
	} , {i,1,(arrayLength)}
								];,
								
	thePoint = Delete[thePoint,0];,
	
	partNumber = Random[Integer,{1,Length[functionList]}],
	functionPart = functionList[[ partNumber ]],
	numberArgs = numList[[ partNumber ]], 
	argumentPart1 = varList[[ Random[Integer,{1,Length[varList]}] ]],
	argumentPart2 = varList[[ Random[Integer,{1,Length[varList]}] ]],
	If[numberArgs == 1,	fullFunction = functionPart[argumentPart1],
				fullFunction = functionPart[argumentPart1,argumentPart2]	
							],
	newIndividual = ReplacePart[exp,fullFunction, {thePoint}]

					 
	}

getPos[exp_] := { Clear[partList];,
	placeCounter=0;, partList[placeCounter] = -1;,
		curExp = exp;, 
		While[Length[curExp] != 0, 
		{++placeCounter,
		partList[placeCounter] = 
					Random[Integer,{1,Length[curExp]}],
		
		curExp = curExp[[partList[placeCounter] ]]
					};],arrayLength = (placeCounter)}

The result of all that is a new array called initPopulation. It has as many members as in popSize (in the function createRandomPop)

The Fitness Test

The one major element that differs from problem to problem is the fitness test—one must be devised for each problem individually. The following function calls fitnessTest for each member of the population and records the score into an array. fitnessTest must produce a global variable called fitnessTestScore, and optionally numHits. numHits is a count of the number succesful attempts. It only applies if the fitnessTest has a timeout or some other type of ceiling.

testAll[popArray_,popSize_,testDepth_, doPrint_]:= {
	rawFitnessArray = {},
	hitsArray = {},
	Do[{Print[popArray[[iter]]],
		fitnessTest[popArray[[iter]],testDepth, False],
		rawFitnessArray = Append[rawFitnessArray, fitnessTestScore],
		hitsArray = Append[hitsArray, numHits],
		If[doPrint == True, {Print["Hits: ", numHits],
							Print[iter,": ", rawFitnessArray[[iter]]	
							 ]} ]}
			,{iter,1,popSize}]
								}

arrayOfOrganisms is the array created by createRandomPop (described above). doPrint is a boolean value which is used to optionally display each score as it is tested. popSize is the size of the population. testAll produces rawFitnessArray which contains the fitness scores of each individual. testDepth is an optional input variable. This is can be used if there are variables in the fitness test itself which may vary form generation to generation.

The following function figures out which function scored the best on the fitness test (best of generation). It also determines the total of scores and the adjusted and normalized fitnesses for each individual . This function assumes that a lower score is better, but it can be easily adjusted if the opposite is true.

analyzeRawFitness[rawFitnessArray_,popSize_] :=	{
	adjustedFitnessArray = {},
	normalizedFitnessArray = {},
	fitnessTotal = rawFitnessArray[[1]],
	bestOfGeneration = rawFitnessArray[[1]],
	bestOfGenNum = 1,
	Do[{fitnessTotal += rawFitnessArray[[i]],
		If[rawFitnessArray[[i]] < bestOfGeneration,{
			bestOfGeneration = rawFitnessArray[[i]],
			bestOfGenNum = i}]}, {i,2,popSize}],
	adjustedTotal = 0,
	Do[{adjustedTotal += (1 / (1 + rawFitnessArray[[i]])),
		adjustedFitnessArray = Append[adjustedFitnessArray, (1 / (1 + rawFitnessArray[[i]]))]},
											{i,1,popSize}],
	Do[normalizedFitnessArray = Append[normalizedFitnessArray, (adjustedFitnessArray[[i]] / adjustedTotal)],
											{i,1,popSize}]
														}

Depending on the type of result rawFitnessArray contains, it may need to be adjusted. Once again, this function assumes lower is better, but it can be easily adjusted (if necessary) by adding a constant to the raw score.

Genetic Crossover

Of course, genetic crossover is the basis of Genetic Algorithms. A bad crossover function will produce bad results for the entire set functions. I have found the following crossover routines themselves to produce very dependable and diverse results for equations with any level of inticracy and size. They have been debugged intensly and revised many times.

The following functions may not be as perfect as the actual crossover routines (which are shown later) but they provide a basis for the batch processing of a generation, after their fitness tests have been completed.

createNewGeneration takes as inputs the last generation computed, that generation's adjusted fitness array, and the number of individuals in the array. All functions necessary for crossover are called from this function.

createNewGeneration[oldGeneration_,oldGenFitArray_,popSize_] := {
numberOfInd = 0,
newGeneration = {},
While[numberOfInd < popSize, {
	pickToCross[oldGenFitArray, popSize],
	If[numberOfPicks == 2, 
		{performCross[oldGeneration[[numberCross1]],oldGeneration[[numberCross2]] ],
		 numberOfInd += 1, newGeneration = Append[newGeneration, newEx1],
		 numberOfInd += 1, newGeneration = Append[newGeneration, newEx2] },
		{numberOfInd += 1, newGeneration = Append[newGeneration, oldGeneration[[numberCross]]]}],
														}]
															}

This next function could probably be improved. I do not have a very good understanding of how new individuals are chosen to breed, so at this point each individual's relative adjusted fitness is simply compared to a random number between 0 and 1 (and if it is lower than the random number it is chosen). There is a one in ten chance that the individual chosen will simply be copied over to the next generation; otherwise it is crossed over with another individual.

pickToCross[arrayOfRelativeFitness_, popSize_] := {
	meth = Random[Integer,{1,10}], 
	If[meth == 1, cross = False, cross = True], 
defined = False,
While[defined == False,{
	rawPick1 = Random[Integer,{1,popSize}],
	If[arrayOfRelativeFitness[[rawPick1]] <= Random[],
		{defined = True, numberCross1 = rawPick1}]
										}],
defined = False,
If[cross == True,
While[defined == False,{
	rawPick2 = Random[Integer,{1,popSize}],
	If[arrayOfRelativeFitness[[rawPick2]] <= Random[],
		{defined = True, numberCross2 = rawPick2}]
										}],
										],
If[cross == True, numberOfPicks = 2, numberOfPicks = 1]

											}

Here are the crossover routines themselves:

performCross[exp1_,exp2_] := {
	If[And[Head[exp1] != Integer,Head[exp2] != Integer,
		   Head[exp1] != Symbol,Head[exp2] != Symbol,
		   Head[exp1] != Complex, Head[exp2] != Complex
		   ],
	{
	doCheck[],
	crossOver[exp1,exp2]
	},
	{newEx1 = exp1, newEx2 = exp2}, {doCheck[], crossOver[exp1,exp2]} ]
	}

The inputs are simply the two equations you want to cross over.
Below are the functions performCross calls. How it works is not important, just realize that it creates two new individuals named newEx1 and newEx2.

doCheck[] := 
	{
	again = True;
	While[again == True, 
		If[Head[exp1] == ReplacePart, exp1 = exp1[[1]], again = False, again = False]
		]
	again = True;
	While[again == True, 
		If[Head[exp2] == ReplacePart, exp2 = exp2[[1]], again = False, again = False] 
		]
	}

crossOver[exp1_,exp2_]:=  {
	Clear[arrayLength,crossOverPointexp2,crossOverPointexp1,point1,point2],
	getCrossPos[exp1];,
	crossOverPointexp1 = {};,
	Do[ { 
	crossOverPointexp1 = Append[crossOverPointexp1,partList[i] ]
	} , {i,1,(arrayLength)}
															];,
		
										
	getCrossPos[exp2];,
	crossOverPointexp2 = {};,
	Do[ { 
	crossOverPointexp2 = Append[crossOverPointexp2,partList[i] ]
	} , {i,1,(arrayLength)}
															];,
												
	crossOverPointexp1 = Delete[crossOverPointexp1,0];,
	crossOverPointexp2 = Delete[crossOverPointexp2,0];,
	
	
	newEx1 = ReplacePart[exp1,exp2 [[crossOverPointexp2]],
					 {crossOverPointexp1}],
	newEx2 = ReplacePart[exp2,exp1 [[crossOverPointexp1]],
					 {crossOverPointexp2}],
															}

getCrossPos[exp_] := { Clear[partList], minimum = 1;,
	placeCounter=0;, partList[placeCounter] = -1;,
		curExp = exp;, 
		While[partList[placeCounter] != 0, 
		{placeCounter = placeCounter + 1,
		
		partList[placeCounter] = 
					Random[Integer,{minimum,Length[curExp]}],
		
		curExp = curExp[[partList[placeCounter] ]],minimum = 0
					};],arrayLength = (placeCounter - 1)};

Example

Cart Centering Problem

The cart centering problem is clearly outlined and explained in [Koza 1]. Briefly, it is an example of a fine control problem, where the computer has complete control of a cart that it can either speed up or slow down. The cart can either go left or right, and the goal is to get it centered on the track and standing still. The inputs to the problem are velocity and position (v and x). The goal is to find the equation that centers the cart in the least amount of time. For this to work properly, the following must be entered:

gt[y_,x_]:= If[y > x, gt[#] = 1, gt[#] = -1];
Unprotect[Abs];
Abs[False] = False;
Abs[True] = True;
Protect[Abs];

The new function gt performs a Lisp-like if function, which gives a result of 1 if true and -1 if false. This is the key equation to be entered as input. Because any function can operate on any argument, Abs[True] and Abs[False] had to be defined. If you think it is more correct, change Abs[False] = False to True. I'm not sure what the absolute value of False is, so I'll leave that up to you to decide.

Here is what initiating the cart problem looks like (but don’t enter it until you have read everything):

funList = {Abs, Divide, Times, Minus, Plus, gt, Power};
argNum =  { 1,	  2,	  2,	 1,		2,	2,	  2};
varList = {x,v};
popSize = 500;
numGens = 50;
levels = 5;
testDepth = 20;

createRandomPop[levels,funList,varList,argNum, popSize, False];

newGeneration = initPopulation;

Do[{
testAll[newGeneration, popSize, testDepth, False],
analyzeRawFitness[rawFitnessArray, popSize],
oldGeneration = newGeneration,
createNewGeneration[oldGeneration, normalizedFitnessArray, popSize],
}, {counter, 1, numGens}]

funList is a list of all the functions to be used in our example, and argNum is a list of the number of arguments required for each function. varlist declares our two main variables. PopSize is the size of the population. This is recomended to be 500 in [Koza 1], but on my computer (Macintosh Quadra 610) the first generation alone took about 36 hours to test. That means about 4.3 minutes for each individual. Keep that in mind when adjusting popSize and numGens (number of generations). (If anyone actually runs this, please let me know how long each individual took to test and what computer you have).

The number of tests for each individual is the number testDepth. After each test is completed, the results are printed. The first number printed (after test #:) is the amount of time it took to center the cart (it will be 10 if the limit was reached). The second number is the position of the cart on the track when it timed out, and the third number is its velocity. If the testDepth = 20 and the popSize = 500, then this will line will be printed 10,000 times per generation. If you actually decide to run this, you may want to comment out the line that prints this to speed up the run (I included this because I often monitor what’s going on.)

Over the course of the tests, it is inevitable that the you will run into an indeterminate expression or two. This fitness test is very good about recognizing bad results and automatically forces a time out when an indeterminate expression comes up.


Here is the actual fitness test for the problem (referenced from testAll):

fitnessTest[thisOrganism_,trialNum_,doPrint_]:= {Clear[testPoint,testResult],
	numHits = 0,
	Do[{testPoint[i] = Random[Integer,{-75,75}],
		testPoint[i] = N[testPoint[i] / 100]
	},
			{i,1,trialNum}],
	Do[{testSpeed[i] = Random[Integer,{-75,75}],
		testSpeed[i] = N[testSpeed[i] / 100]
	},
			{i,1,trialNum}],
			
	Do[{thePoint = testPoint[i], theDirection = testSpeed[i],
		time = 0,
		thePlace = thePoint,
		theSpeed = theDirection,
	While[((thePlace != 0.0) || (theSpeed != 0.0)) && (time < 10),
		{Clear[x,v],
		time += 0.02,
		modifyDirection = thisOrganism /.{x->thePoint,v->theDirection},
		If[modifyDirection > 0, modifyDirection = 1, modifyDirection = -1,time = 10],
		theDirection += (modifyDirection * 0.02),
		theDirection = Chop[theDirection],
		thePoint += (theDirection * 0.02), 
		thePoint = Chop[thePoint],
		If[doPrint == True, 
			Print[time, " ",thePoint," ", theDirection," ", modifyDirection] ],
		thePlace = thePoint,
		If[thePlace <= 0.01 && thePlace >= -0.01, thePlace = 0.0, {}, time = 10],
		theSpeed = theDirection,
		If[theSpeed < 0.015 && theSpeed > -0.015, theSpeed = 0.0, {}, time = 10]
							}],
				testResult[i] = time,
				If[time != 10., numHits += 1],
				Print["Test ",i,": ",testResult[i], " ",thePoint, " ", theDirection]},	
								{i,1,trialNum}],
				
		totalTime = testResult[1];,
		Do[totalTime += testResult[i], {i,2,trialNum}],
		fitnessTestScore = totalTime / trialNum;,
		Print[fitnessTestScore]
											}

This is not a perfect, but it provides results very similar to those in [Koza 1]. I wasn't sure how much leeway I could give the cart to be considered perfectly centered, so I decided to use the two numbers that gave the closest results to the book.

Contacting the Author

If you have any questions, comments, or suggestions, please e-mail me at: SonOfaBitS@aol.com. I would be very happy to clarify or discuss anything in this document.

If you want to distribute this outside of MathSource, please contact me first.

Bibliography

Koza 1:
Koza, John.1992. Genetic Programming: on the programming of computers by means of natural selection. MIT Press.

Gray:
Gray, John. 1994. Mastering Mathematica.. Academic Press, Inc.

Abelson:
Abelson, H. and Sussman, G.J. with Sussman, J. 1985. Structure and Interpretation of Computer Programs. The MIT Press.