From: schumach@sol.UVic.CA (Ian Schumacher)
Subject: genetic algorithm example in c++
Organization: University of Victoria
Distribution: all
Date: 31 Aug 92 00:02:19 GMT


I have seen a lot of interest in Gene algorithms these days.  I was
pretty interested myself, so I threw together a Gene and Gene_Pool
Class in C++.  Figured I'd send it out, so that anyone that was interested
could have a example to work with, improve upon, whatever.
Bugs?  You tell me.  It hasn't been very thoroughly tested, but there's
not that much that could go wrong.  It's in three parts gene.h gene.m and
imp.m.  Enjoy.

--------------------------------------------------------------------

// class gene.h: gene structure
// Programmer : Ian Schumacher
// Comments welcome
// email: schumach@sol.UVic.CA
// Genie : I.SCHUMACHER

class Gene
{
	unsigned short *gene;

protected:

	int length;
	double mutate_rate;
	double value;

public:
	Gene(int);
	Gene();
	~Gene();
	Gene& operator+(Gene &);
	unsigned short& operator[](int);
	Gene& operator=(const Gene &);
	Gene& operator=(const double);
	Gene& mutate(void);
	void set_rate(double);
	int get_length(void);
	double get_value(void) {return value;}
};

class Gene_Pool
{
	int size;
	int plength;
	Gene *gene;
	int *index;
	
public:
	Gene_Pool(int number=1,int a=1);
	~Gene_Pool() {delete gene;delete index;}
	Gene_Pool& operator+(Gene_Pool &);
	Gene& operator[](int);
	Gene_Pool& operator=(const Gene_Pool &);
	Gene_Pool& mutate(void);
	int get_size(void) { return size;}
	int get_length(void) {return plength;}
	Gene_Pool& evolve(void);
	Gene_Pool& sort(void);
	Gene_Pool& kill(int);
	Gene_Pool& set_rate(double);
};

---------------------------------------------------------------------------

// gene.m : Definition of Gene and Gene_Pool functions
// Programmer : Ian Schumacher
// Comments welcome
// email: schumach@sol.UVic.CA
// Genie : I.SCHUMACHER

#include "gene.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define RATE 0.01			// default mutation rate
#define RAND_MATE 1		// mate randomly or in order switch

void indexx(int ,double *,int *);
int irbit2(void);

int pop_length=1;		// gene population value set by Gene_Pool, used by Gene
int iseed=-666;               // random generator seed

Gene::Gene(int n)
{
	int i;
	
	if(n>0)
	{
		gene=new unsigned short[n];
		mutate_rate=RATE;
		length=n;
	}
	else
	{
		gene=new unsigned short[1];
		mutate_rate=RATE;
		length=1;
	}
	
// insert random values into gene
	for(i=0;i<length;i++)
		gene[i]=(unsigned short) irbit2();  // don't trust rand() that much, use "NR in C" function
}

Gene::Gene()
{
// uses global variable pop_length set by Gene_Pool to determine gene size
	int i;
	
	if(pop_length<1)
		pop_length=1;
	length=pop_length;
	mutate_rate=0.01;
	gene=new unsigned short[pop_length];
// insert random values into gene
	for(i=0;i<pop_length;i++)
		gene[i]=(unsigned short) (irbit2());
}

Gene::~Gene()
{
	delete gene;
}

Gene& Gene::operator+(Gene &a)
{
	int cp,i;
	int tmp;

	static Gene b(length);	//static because its adress is returned

	if(this==&a)
		return *this;
	if(length!=a.length)
	{
		printf("Currently cannot mate genes of different lengths. Genes unchanged\n");
		return *this;
	}
	// mating routine
	cp=int(rand()%length);   // Random cross point. I thought this would produce better result
						// constant half way cross point.
	for(i=0;i<cp;i++)
		b.gene[i]=this->gene[i];
	for(i=cp;i<length;i++)
		b.gene[i]=a.gene[i];
	return b;
}

unsigned short& Gene::operator[](int elem)
// Returning address allows gene[j] to be on either side of equals sign
{
	static unsigned short tmp=0;
	
	if(elem>=0&&elem<length)
		return gene[elem];
	else
	{
		printf("Gene access error\n");
		return tmp;
	}
}

Gene& Gene::operator=(const Gene& a)
// set one gene equal to another
{
	int i;

	if(this==&a)
		return *this;
	length=a.length;
	mutate_rate=a.mutate_rate;
	value=a.value;
	delete gene;
	gene=new unsigned short[length];
	for(i=0;i<length;i++)
		gene[i]=a.gene[i];
	return *this;
}

Gene& Gene::operator=(const double x)
// set gene value
{
	value=x;
	return *this;
}

Gene& Gene::mutate(void)
{
	int i;
	
	for(i=0;i<length;i++)
	{
		if(rand()/(RAND_MAX+1.0)<mutate_rate)
		{
			gene[i]=gene[i]<1 ? 1:0;
		}
	}
	return *this;
}

void Gene::set_rate(double r)
{
	if(r>0.0)
		mutate_rate=r;
	return;
}

int Gene::get_length(void)
{
	return length;
}

// Gene_Pool definitions

Gene_Pool::Gene_Pool(int number,int a)
// uses global variable pop_length to set gene size
{
	int i;
	
	pop_length=a;
	plength=a;
	size=number;
	index=new int[size];
	for(i=0;i<size;i++)
		index[i]=i;
	gene=new Gene[size];
}

Gene_Pool& Gene_Pool::operator+(Gene_Pool &a)
{
	int i;
	
	static Gene_Pool b(size+a.size,plength);
	
	if(this==&a)
	{
		delete &b;
		return *this;
	}
	if(plength!=a.plength)
	{
		printf("incompatible genes\n");
		delete &b;
		return *this;
	}
	for(i=0;i<a.size;i++)
		b.gene[i]=gene[i];
	for(i=size;i<b.size;i++)
		b.gene[i]=a.gene[i-size];
	return b;
}

Gene& Gene_Pool::operator[](int elem)
{
	if(elem>=0&&elem<size)
		return gene[index[size-elem-1]];
	else
	{
		printf("Gene access error\n");
		return gene[0];
	}
}

Gene_Pool& Gene_Pool::operator=(const Gene_Pool& a)
{
	int i;

	if(this==&a)
		return *this;
	size=a.size;
	pop_length=a.plength;
	plength=a.plength;
	delete gene;
	delete index;
	gene=new Gene[size];
	index=new int[size];
	for(i=0;i<pop_length;i++)
	{
		gene[i]=a.gene[i];
		index[i]=a.index[i];
	}
	return *this;
}

Gene_Pool& Gene_Pool::mutate(void)
{
	int i;

	for(i=0;i<size;i++)
		gene[i].mutate();
	return *this;
}

Gene_Pool& Gene_Pool::set_rate(double r)
{
	int i;

	for(i=0;i<size;i++)
		gene[i].set_rate(r);
	return *this;
}

Gene_Pool& Gene_Pool::sort(void)
{
	int i,j;
	double *values;

	values=new double[size];

	for(i=0;i<size;i++)
		values[i]=gene[i].get_value();
	indexx(size,values,index);	// "NR in C" heap sort function.  Sorts lowest to highest
	delete values;
	return *this;
}

Gene_Pool& Gene_Pool::evolve(void)
// allows for random mating of top half or ordered mating (you decide which is better)
{
	int i;

	this->sort();
	if(RAND_MATE)
	{
		// randomly mate top half of population and store in bottom half
		for(i=0;i<size/2;i++)
		{
			gene[index[i]]=gene[index[size-rand()%size-1]]+gene[index[size-rand()%size-1]];
		}
		this->mutate();
	}
	else
	{
		// or mate top half of population together in order
		for(i=0;i<size/2;i++)
		{
			gene[index[i]]=gene[index[size-i-1]]+gene[index[size-i-2]];
		}
	}
	return *this;
}

Gene_Pool& Gene_Pool::kill(int n)
// if after evolve for a while you want to kill off the weak to speed things up
{
	int i;
	Gene *tmp;

	tmp=new Gene[size-n];
	if(n>=size)
		n=size-1;
	size=size-n;
	for(i=0;i<size;i++)
		tmp[i]=gene[index[size+n-i]-1];
	delete gene;
	delete index;
	gene=new Gene[size];
	index=new int[size];
	for(i=0;i<size;i++)
		gene[i]=tmp[i];
	this->sort();
	delete tmp;
	return *this;
}

// "Numerical Recipes is C" heapsort of arrin with corresponding reordering of indx
void indexx(int n,double *arrin,int *indx)
{
	int l,j,ir,indxt,i;
	float q;

	indx--;
	for (j=1;j<=n;j++) indx[j]=j-1;
	l=(n >> 1) + 1;
	ir=n;
	for (;;)
	{
		if (l > 1)
			q=arrin[(indxt=indx[--l])];
		else
		{
			q=arrin[(indxt=indx[ir])];
			indx[ir]=indx[1];
			if (--ir == 1)
			{
				indx[1]=indxt;
				indx++;
				return;
			}
		}
		i=l;
		j=l << 1;
		while (j <= ir)
		{
			if (j < ir && arrin[indx[j]] < arrin[indx[j+1]]) j++;
			if (q < arrin[indx[j]])
			{
				indx[i]=indx[j];
				j += (i=j);
			}
			else j=ir+1;
		}
		indx[i]=indxt;
	}
}

#define IB1 1
#define IB2 2
#define IB5 16
#define IB18 131072
#define MASK IB1+IB2+IB5

// "Numerical Recipes in C" random bit function (uses periodic maximal length sequences)
int irbit2(void)
{
	if (
	iseed & IB18) {
		iseed=((iseed ^ MASK) << 1) | IB1;
		return 1;
	} else {
		iseed <<= 1;
		return 0;
	}
}

#undef MASK
#undef IB18
#undef IB5
#undef IB2
#undef IB1

--------------------------------------------------------------------------

// imp.m : Example implementation of Gene_Pool class
// Programmer : Ian Schumacher
// Comments welcome 
// email : schumach@sol.UVic.CA
// Genie : I.SCHUMACHER

// Pathetically unrealistic problem to be solved by using gene algorithms,
// but illustrates principles.

#include "gene.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define size 500		// Gene_Pool population size
#define length 30		// Gene length
#define steps 100		// maximum number of evolution steps
#define sqr(x) x*x		// simple squaring function

int main()
{
	int i,j,k,t;
	Gene_Pool a(size,length);
	double value1,value2,value3;
	double tmp1,tmp2;
	int tmp;

	a.set_rate(0.05/length);		// set mutation rate, (in this case to 5% of population)

	for(i=0;i<steps;i++)		// number of evolution steps
	{
		for(j=0;j<a.get_size();j++)
		{
			// this is where the value of gene is determined
			value1=0;
			value2=0;
			value3=0;
			tmp=a[j].get_length()/3;
			for(t=0;t<tmp;t++)
			{
				// interperate gene as three long integers
				value1+=(1L<<t)*a[j][t];
				value2+=(1L<<t)*a[j][t+tmp];
				value3+=(1L<<t)*a[j][t+2*tmp];
			}
			// normalize values to 1
			value1/=1L<<tmp;
			value2/=1L<<tmp;
			value3/=1L<<tmp;
			// maximize simple nonlinear equation
			a[j]=exp(-sqr((value1-0.33))-sqr((value2-0.5))-sqr((value3-0.7)));  // set gene[j] value
		}
// break out of evolution early if genes starting to look the same
		if((a[0].get_value()==a[1].get_value())&&a[1].get_value()==a[2].get_value()&&tmp1==tmp2)
			break;
		tmp1=a[0].get_value();
		tmp2=a[1].get_value();
		a.evolve();
		// display sample of every fifth generation
		if(!(i%5))
		{
			for(j=0;j<5;j++)
			{
				for(k=0;k<length;k++)
				{
					printf("%d ",a[j][k]);
				}
				printf("%f\n",a[j].get_value());
			}
			printf("\n");
		}
	}
	// display final result
	for(t=0;t<tmp;t++)
	{
		value1+=(1L<<t)*a[0][t];
		value2+=(1L<<t)*a[0][t+tmp];
		value3+=(1L<<t)*a[0][t+2*tmp];
	}
	value1/=1L<<tmp;
	value2/=1L<<tmp;
	value3/=1L<<tmp;
	printf("value1 = %f value2= %f value3= %f\n",value1,value2,value3);
	return 1;
}