Differential Evolution C++ Code - University of California, Berkeley

DESolver.cpp

#include <memory.h>
#include "DESolver.h"

#define Element(a,b,c)  a[b*nDim+c]
#define RowVector(a,b)  (&a[b*nDim])
#define CopyVector(a,b) memcpy((a),(b),nDim*sizeof(double))

DESolver::DESolver(int dim,int popSize) :
					nDim(dim), nPop(popSize),
					generations(0), strategy(stRand1Exp),
					scale(0.7), probability(0.5), bestEnergy(0.0),
					trialSolution(0), bestSolution(0),
					popEnergy(0), population(0)
{
	trialSolution = new double[nDim];
	bestSolution  = new double[nDim];
	popEnergy	  = new double[nPop];
	population	  = new double[nPop * nDim];

	return;
}

DESolver::~DESolver(void)
{
	if (trialSolution) delete trialSolution;
	if (bestSolution) delete bestSolution;
	if (popEnergy) delete popEnergy;
	if (population) delete population;

	trialSolution = bestSolution = popEnergy = population = 0;
	return;
}

void DESolver::Setup(double *min,double *max,
						int deStrategy,double diffScale,double crossoverProb)
{
	int i;

	strategy	= deStrategy;
	scale		= diffScale;
	probability = crossoverProb;
	
	for (i=0; i < nPop; i++)
	{
		for (int j=0; j < nDim; j++)
			Element(population,i,j) = RandomUniform(min[j],max[j]);

		popEnergy[i] = 1.0E20;
	}

	for (i=0; i < nDim; i++)
		bestSolution[i] = 0.0;

	switch (strategy)
	{
		case stBest1Exp:
			calcTrialSolution = Best1Exp;
			break;

		case stRand1Exp:
			calcTrialSolution = Rand1Exp;
			break;

		case stRandToBest1Exp:
			calcTrialSolution = RandToBest1Exp;
			break;

		case stBest2Exp:
			calcTrialSolution = Best2Exp;
			break;

		case stRand2Exp:
			calcTrialSolution = Rand2Exp;
			break;

		case stBest1Bin:
			calcTrialSolution = Best1Bin;
			break;

		case stRand1Bin:
			calcTrialSolution = Rand1Bin;
			break;

		case stRandToBest1Bin:
			calcTrialSolution = RandToBest1Bin;
			break;

		case stBest2Bin:
			calcTrialSolution = Best2Bin;
			break;

		case stRand2Bin:
			calcTrialSolution = Rand2Bin;
			break;
	}

	return;
}

bool DESolver::Solve(int maxGenerations)
{
	int generation;
	int candidate;
	bool bAtSolution;

	bestEnergy = 1.0E20;
	bAtSolution = false;

	for (generation=0;(generation < maxGenerations) && !bAtSolution;generation++)
		for (candidate=0; candidate < nPop; candidate++)
		{
			(this->*calcTrialSolution)(candidate);
			trialEnergy = EnergyFunction(trialSolution,bAtSolution);

			if (trialEnergy < popEnergy[candidate])
			{
				// New low for this candidate
				popEnergy[candidate] = trialEnergy;
				CopyVector(RowVector(population,candidate),trialSolution);

				// Check if all-time low
				if (trialEnergy < bestEnergy)
				{
					bestEnergy = trialEnergy;
					CopyVector(bestSolution,trialSolution);
				}
			}
		}

	generations = generation;
	return(bAtSolution);
}

void DESolver::SelectSamples(int candidate,int *r1,int *r2,
										int *r3,int *r4,int *r5)
{
	if (r1)
	{
		do
		{
			*r1 = (int)RandomUniform(0.0,(double)nPop);
		}
		while (*r1 == candidate);
	}

	if (r2)
	{
		do
		{
			*r2 = (int)RandomUniform(0.0,(double)nPop);
		}
		while ((*r2 == candidate) || (*r2 == *r1));
	}

	if (r3)
	{
		do
		{
			*r3 = (int)RandomUniform(0.0,(double)nPop);
		}
		while ((*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1));
	}

	if (r4)
	{
		do
		{
			*r4 = (int)RandomUniform(0.0,(double)nPop);
		}
		while ((*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1));
	}

	if (r5)
	{
		do
		{
			*r5 = (int)RandomUniform(0.0,(double)nPop);
		}
		while ((*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3)
													|| (*r5 == *r2) || (*r5 == *r1));
	}

	return;
}

/*------------  Strategies  ---------------------------------------------*/

void DESolver::Best1Exp(int candidate)
{
	int r1, r2;
	int n;

	SelectSamples(candidate,&r1,&r2);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
	{
		trialSolution[n] = bestSolution[n]
							+ scale * (Element(population,r1,n)
							- Element(population,r2,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Rand1Exp(int candidate)
{
	int r1, r2, r3;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
	{
		trialSolution[n] = Element(population,r1,n)
							+ scale * (Element(population,r2,n)
							- Element(population,r3,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::RandToBest1Exp(int candidate)
{
	int r1, r2;
	int n;

	SelectSamples(candidate,&r1,&r2);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
	{
		trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
							 + scale * (Element(population,r1,n)
							 - Element(population,r2,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Best2Exp(int candidate)
{
	int r1, r2, r3, r4;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3,&r4);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
	{
		trialSolution[n] = bestSolution[n] +
							scale * (Element(population,r1,n)
										+ Element(population,r2,n)
										- Element(population,r3,n)
										- Element(population,r4,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Rand2Exp(int candidate)
{
	int r1, r2, r3, r4, r5;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3,&r4,&r5);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
	{
		trialSolution[n] = Element(population,r1,n)
							+ scale * (Element(population,r2,n)
										+ Element(population,r3,n)
										- Element(population,r4,n)
										- Element(population,r5,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Best1Bin(int candidate)
{
	int r1, r2;
	int n;

	SelectSamples(candidate,&r1,&r2);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; i < nDim; i++) 
	{
		if ((RandomUniform(0.0,1.0) < probability) || (i == (nDim - 1)))
			trialSolution[n] = bestSolution[n]
								+ scale * (Element(population,r1,n)
											- Element(population,r2,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Rand1Bin(int candidate)
{
	int r1, r2, r3;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; i < nDim; i++) 
	{
		if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
			trialSolution[n] = Element(population,r1,n)
								+ scale * (Element(population,r2,n)
												- Element(population,r3,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::RandToBest1Bin(int candidate)
{
	int r1, r2;
	int n;

	SelectSamples(candidate,&r1,&r2);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; i < nDim; i++) 
	{
		if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
			trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
									+ scale * (Element(population,r1,n)
												- Element(population,r2,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Best2Bin(int candidate)
{
	int r1, r2, r3, r4;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3,&r4);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; i < nDim; i++) 
	{
		if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
			trialSolution[n] = bestSolution[n]
								+ scale * (Element(population,r1,n)
											+ Element(population,r2,n)
											- Element(population,r3,n)
											- Element(population,r4,n));
		n = (n + 1) % nDim;
	}

	return;
}

void DESolver::Rand2Bin(int candidate)
{
	int r1, r2, r3, r4, r5;
	int n;

	SelectSamples(candidate,&r1,&r2,&r3,&r4,&r5);
	n = (int)RandomUniform(0.0,(double)nDim);

	CopyVector(trialSolution,RowVector(population,candidate));
	for (int i=0; i < nDim; i++) 
	{
		if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
			trialSolution[n] = Element(population,r1,n)
								+ scale * (Element(population,r2,n)
											+ Element(population,r3,n)
											- Element(population,r4,n)
											- Element(population,r5,n));
		n = (n + 1) % nDim;
	}

	return;
}

/*------Constants for RandomUniform()---------------------------------------*/
#define SEED 3
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double DESolver::RandomUniform(double minValue,double maxValue)
{
	long j;
	long k;
	static long idum;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	double result;

	if (iy == 0)
		idum = SEED;

	if (idum <= 0)
	{
		if (-idum < 1)
			idum = 1;
		else
			idum = -idum;

		idum2 = idum;

		for (j=NTAB+7; j>=0; j--)
		{
			k = idum / IQ1;
			idum = IA1 * (idum - k*IQ1) - k*IR1;
			if (idum < 0) idum += IM1;
			if (j < NTAB) iv[j] = idum;
		}

		iy = iv[0];
	}

	k = idum / IQ1;
	idum = IA1 * (idum - k*IQ1) - k*IR1;

	if (idum < 0)
		idum += IM1;

	k = idum2 / IQ2;
	idum2 = IA2 * (idum2 - k*IQ2) - k*IR2;

	if (idum2 < 0)
		idum2 += IM2;

	j = iy / NDIV;
	iy = iv[j] - idum2;
	iv[j] = idum;

	if (iy < 1)
		iy += IMM1;

	result = AM * iy;

	if (result > RNMX)
		result = RNMX;

	result = minValue + result * (maxValue - minValue);
	return(result);
}

DESolver.h

#if !defined(_DESOLVER_H)
#define _DESOLVER_H

#define stBest1Exp			0
#define stRand1Exp			1
#define stRandToBest1Exp	2
#define stBest2Exp			3
#define stRand2Exp			4
#define stBest1Bin			5
#define stRand1Bin			6
#define stRandToBest1Bin	7
#define stBest2Bin			8
#define stRand2Bin			9

class DESolver;

typedef void (DESolver::*StrategyFunction)(int);

class DESolver
{
public:
	DESolver(int dim,int popSize);
	~DESolver(void);
	
	// Setup() must be called before solve to set min, max, strategy etc.
	void Setup(double min[],double max[],int deStrategy,
							double diffScale,double crossoverProb);

	// Solve() returns true if EnergyFunction() returns true.
	// Otherwise it runs maxGenerations generations and returns false.
	virtual bool Solve(int maxGenerations);

	// EnergyFunction must be overridden for problem to solve
	// testSolution[] is nDim array for a candidate solution
	// setting bAtSolution = true indicates solution is found
	// and Solve() immediately returns true.
	virtual double EnergyFunction(double testSolution[],bool &bAtSolution) = 0;
	
	int Dimension(void) { return(nDim); }
	int Population(void) { return(nPop); }

	// Call these functions after Solve() to get results.
	double Energy(void) { return(bestEnergy); }
	double *Solution(void) { return(bestSolution); }

	int Generations(void) { return(generations); }

protected:
	void SelectSamples(int candidate,int *r1,int *r2=0,int *r3=0,
												int *r4=0,int *r5=0);
	double RandomUniform(double min,double max);

	int nDim;
	int nPop;
	int generations;

	int strategy;
	StrategyFunction calcTrialSolution;
	double scale;
	double probability;

	double trialEnergy;
	double bestEnergy;

	double *trialSolution;
	double *bestSolution;
	double *popEnergy;
	double *population;

private:
	void Best1Exp(int candidate);
	void Rand1Exp(int candidate);
	void RandToBest1Exp(int candidate);
	void Best2Exp(int candidate);
	void Rand2Exp(int candidate);
	void Best1Bin(int candidate);
	void Rand1Bin(int candidate);
	void RandToBest1Bin(int candidate);
	void Best2Bin(int candidate);
	void Rand2Bin(int candidate);
};

#endif // _DESOLVER_H

DETest.cpp

#include <stdio.h>
#include "DESolver.h"

// Polynomial fitting problem
class PolynomialSolver : public DESolver
{
public:
	PolynomialSolver(int dim,int pop) : DESolver(dim,pop), count(0) {;}
	double EnergyFunction(double trial[],bool &bAtSolution);

private:
	int count;
};

double PolynomialSolver::EnergyFunction(double *trial,bool &bAtSolution)
{
	int i, j;
	int const M=60;
	double px, x=-1, dx=M, result=0;

	dx = 2.0 / dx;
	for (i=0; i<=M; i++)
	{
		px = trial[0];
		for (j=1;j<nDim;j++)
			px = x*px + trial[j];

		if (px<-1 || px>1)
			result += (1 - px) * (1 - px);

		x += dx;
	}

	px = trial[0];
	for (j=1;j<nDim;j++)
		px = 1.2*px + trial[j];

	px = px - 72.661;
	if (px<0)
		result += px * px;

	px = trial[0];
	for (j=1; j<nDim; j++)
		px = -1.2*px + trial[j];

	px = px - 72.661;
	
	if (px<0)
		result+=px*px;

	if (count++ % nPop == 0)
		printf("%d %lf\n",count / nPop + 1,Energy());
	
	return(result);
}

#define N_DIM 9
#define N_POP 100
#define MAX_GENERATIONS	800
int main(void)
{
	double min[N_DIM];
	double max[N_DIM];
	int i;

	PolynomialSolver solver(N_DIM,N_POP);

	for (i=0;i<N_DIM;i++)
	{
		max[i] =  100.0;
		min[i] = -100.0;
	}

	solver.Setup(min,max,stBest1Exp,0.9,1.0);
	
	printf("Calculating...\n\n");
	solver.Solve(MAX_GENERATIONS);

	double *solution = solver.Solution();

	printf("\n\nBest Coefficients:\n");
	for (i=0;i<N_DIM;i++)
		printf("[%d]: %lf\n",i,solution[i]);

	return;
}