// FILENAME: EvlStrat.cpp
// AUTHOR: Andrew Volz. Note that the main evolution strategy algorithm, found in
//	methods generateNewIndividual and evolve, have been adapted from code supplied
//	in the book "Evolution and Optimum Seeking" (1995) by Hans-Paul Schwefel.
// CREATED: 23/08/96
// LAST MODIFIED: 24/08/96
// DESCRIPTION: Implementation of the EvolutionStrategy class.



// Include files
#include "stdafx.h"
#include "EvlStrat.h"
#include <stdlib.h>
#include <math.h>
#include <time.h>



// FUNCTION: Default Constructor
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 23/08/96
// DESCRIPTION: Default constructor for EvolutionStrategy class.

EvolutionStrategy::EvolutionStrategy()
{
	noOfIterations = 0;
	noOfMutations = 0;
	srand(time(NULL));
	currentParentFitness = 0;
	currentOffspringFitness = 0;
}



// FUNCTION: Constructor
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 23/08/96
// DESCRIPTION: Constructor for EvolutionStrategy class.  This constructor should
//	be used when the parameter settings of the evolution strategy class are known
//	before the EvolutionStrategy class is created.

EvolutionStrategy::EvolutionStrategy(const ESConfiguration &intendedConfiguration)
{
	noOfIterations = 0;
	noOfMutations = 0;
	configureParameters(intendedConfiguration);
	copyConfiguration(initialConfiguration);
	srand(time(NULL));
	currentParentFitness = 0;
	currentOffspringFitness = 0;
}



// FUNCTION: Destructor
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 23/08/96
// DESCRIPTION: Destructor for EvolutionStrategy class.

EvolutionStrategy::~EvolutionStrategy()
{
}



// FUNCTION: configureParameters
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 23/08/96
// DESCRIPTION: Configures the evolution strategy parameters based on
//		the intendedConfiguration structure.  If a parameter is invalid,
//		it is set to it's default value.  Parameters xb and sm must not be
//		empty.

int EvolutionStrategy::configureParameters (const ESConfiguration &intendedConfiguration)
{
	int parametersOK = TRUE;

	if (intendedConfiguration.n > minNValue)
		initialConfiguration.n = intendedConfiguration.n;
	else
		initialConfiguration.n = nDefault;
	if (intendedConfiguration.m >= minMValue)
		initialConfiguration.m = intendedConfiguration.m;
	else
		initialConfiguration.m = mDefault;
	initialConfiguration.lf = 0;
	initialConfiguration.lr = intendedConfiguration.lr;
	if (intendedConfiguration.ls >= minLsValue)
		initialConfiguration.ls = intendedConfiguration.ls;
	else
		initialConfiguration.ls = lsDefault;
	if (intendedConfiguration.tm > minTmValue)
		initialConfiguration.tm = intendedConfiguration.tm;
	else
		initialConfiguration.tm = tmDefault;
	initialConfiguration.ea = intendedConfiguration.ea;
	initialConfiguration.eb = intendedConfiguration.eb;
	initialConfiguration.ec = intendedConfiguration.ec;
	initialConfiguration.ed = intendedConfiguration.ed;
	initialConfiguration.sn = intendedConfiguration.sn;
	initialConfiguration.fb = 0.0;
	copyParameters(initialConfiguration.xb, intendedConfiguration.xb, intendedConfiguration.n+1);
	copyParameters(initialConfiguration.sm, intendedConfiguration.sm, intendedConfiguration.n+1);
	copyConfiguration(initialConfiguration);

	return(parametersOK);
}



// FUNCTION: copyParameters
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 26/08/96
// DESCRIPTION: Copies parameterList2 to parameterList1, up to the
//		ith subscript.

void EvolutionStrategy::copyParameters(double *parameterList1, const double *parameterList2, const unsigned int &i)
{
	for (unsigned int counter=0; counter < i; counter++)
		parameterList1[counter] = parameterList2[counter];
	return;
}



// FUNCTION: initialiseStrategy
// AUTHOR: Adapted from Schwefel's (1995) _evol function by Andrew Volz
// CREATED: 24/08/96
// LAST MODIFIED: 24/08/96
// DESCRIPTION: Initialises the evolution strategy parameters which are
//		NOT provided by the user. This must be run before the
//		generateNewPopulation function may be executed (but this
//		is handled by the EvolutionStrategy constructor).

int EvolutionStrategy::initialiseStrategy()
{
	evz_1.lz = 1L;
	if (currentConfiguration.m <= 0)
	{
		// Problem has no restrictions
		currentConfiguration.lf = 1L;
		currentConfiguration.fb = objectiveFunction(currentConfiguration.xb);
		currentParentFitness = currentConfiguration.fb;
		for (k = 1; k <= 10; ++k)
			l[k - 1] = currentConfiguration.n * k / 5L;
		le = currentConfiguration.n + currentConfiguration.n;
		lm = 0L;
		lc = 0L;
		fc = currentConfiguration.fb;
	}
	else
	{
		// Handle restrictions
		currentConfiguration.lf = -1L;
		currentConfiguration.fb = (double) 0.;
		i_1 = currentConfiguration.m;
		for (j = 1; j <= i_1; ++j)
		{
			fg = constraintsFunction(j, currentConfiguration.n, currentConfiguration.xb);
			if (fg < 0.)
				currentConfiguration.fb -= fg;	// Decrement current best
		}
		if (currentConfiguration.fb <= 0.)
		{
			currentConfiguration.lf = 1L;
			currentConfiguration.fb = objectiveFunction(currentConfiguration.xb);
			currentParentFitness = currentConfiguration.fb;
		}
		for (k = 1; k <= 10; ++k)
			l[k - 1] = currentConfiguration.n * k / 5L;
		le = currentConfiguration.n + currentConfiguration.n;
		lm = 0L;
		lc = 0L;
		fc = currentConfiguration.fb;
	}
	return(currentStatus=READY);
}



// FUNCTION: generateNewIndividual
// AUTHOR: Adapted from Schwefel's (1995) _evol function by Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 24/08/96
// DESCRIPTION: Executes a single iteration of the evolution strategy
//		process, ie. generates one new population of solutions. Returns
//		error code.

long int EvolutionStrategy::generateNewIndividual()
{
	for (int counter = 1; counter <= currentConfiguration.n; counter++)
		currentConfiguration.x[counter] = currentConfiguration.xb[counter] + gaussianRandomNumber(currentConfiguration.sm[counter]);
	noOfMutations++;
	return(TRUE);
}



// FUNCTION: gaussianRandomNumber
// AUTHOR: Modified from Schwefel's code by Andrew Volz
// CREATED: Unknown
// LAST MODIFIED: 26/08/96
// DESCRIPTION: Returns a Gaussian random number using Box-Mueller transformation rules.

double EvolutionStrategy::gaussianRandomNumber(double s)
{
		static double zp = 6.28318531;
		double ret_val;
		static double a, b, d;

		switch((int) evz_1.lz)
		{
		case 1:
			a = sqrt(log(uniformRandomNumber()) * (double) -2.);
			b = zp * uniformRandomNumber();
			ret_val = s * a * sin(b);
			evz_1.lz = 2L;
			break;
		default:
			ret_val = s *a * cos(b);
			evz_1.lz = 1L;
			break;
		}

		return(ret_val);
}



// FUNCTION: resetStrategy
// AUTHOR: Andrew Volz
// CREATED: 23/08/96
// LAST MODIFIED: 24/08/96
// DESCRIPTION: Resests all evolution strategy parameters to their default values.

int EvolutionStrategy::resetStrategy()
{
	currentParentFitness = 0;
	currentOffspringFitness = 0;
	noOfIterations = 0;
	noOfMutations = 0;
	return(TRUE);
}



// FUNCTION: iterate
// AUTHOR: Andrew Volz
// CREATED: 24/08/96
// LAST MODIFIED: 06/09/96
// DESCRIPTION: Either generates a new individual, or chooses the best of the two
//		individuals and then generates offspring from the best, depending
//		on the stage of the evolution strategy process.  Returns an array
//		containing the modified parameters.

int EvolutionStrategy::iterate(double xVector[])
{
	if(noOfIterations==0)
		initialiseStrategy();

	evolve();
	if (currentStatus == INCOMPLETE)
		copyParameters(xVector, currentConfiguration.x, currentConfiguration.n+1);
	else
		copyParameters(xVector, currentConfiguration.xb, currentConfiguration.n+1);

	noOfIterations++;
	return(currentStatus);
}



// FUNCTION: evolve
// AUTHOR: Adapted from Shwefel's _evol function by Andrew Volz
// CREATED: 25/08/96
// LAST MODIFIED: 06/09/96
// DESCRIPTION: Chooses current best solution.  Generates offspring from this solution,
//		(offspring held in currentConfiguration.x).

int EvolutionStrategy::evolve()
{
		if (currentStatus == INCOMPLETE)
			goto L15;

	L7:	
		generateNewIndividual();
	
		if (currentConfiguration.lf <= 0)
		{
			// Starting point is not feasible
			ff = (double) 0.;
			i_1 = currentConfiguration.m;
			for (j = 1; j <= i_1; ++j)
			{
				fg = constraintsFunction(j, currentConfiguration.n, currentConfiguration.x);
				if (fg < 0.)
					ff -= fg;
			}
			if (ff <= 0.)
			{
				// All constraints have been satisfied
				i_1 = currentConfiguration.n;
				for (i = 1; i <= i_1; ++i)
					currentConfiguration.xb[i] = currentConfiguration.x[i];
				currentConfiguration.fb = objectiveFunction(currentConfiguration.xb);
				currentParentFitness = currentConfiguration.fb;
				currentConfiguration.lf = 0L;
				// Starting point not feasible, search for a feasible state
				// succesful.
				return(currentStatus=COMPLETE);
			}
			else
				goto L16;
		}

		
		if (currentConfiguration.m > 0)
		{
			// Handle constraints
			i_1 = currentConfiguration.m;
			for (j = 1; j <= i_1; ++j)
			{
				if (constraintsFunction(j, currentConfiguration.n, currentConfiguration.x) < 0.)
					goto L19;
			}
		}

		return(currentStatus=INCOMPLETE);

	L15:
		ff = objectiveFunction(currentConfiguration.x);
		currentOffspringFitness = ff;

	L16:
		if (ff - currentConfiguration.fb <= 0.)
		{
			// Offspring is fitter than parent, update xb appropriately.
			++le;
			currentConfiguration.fb = ff;
			i_1 = currentConfiguration.n;
			for (i = 1; i <= i_1; ++i)
				currentConfiguration.xb[i] = currentConfiguration.x[i];
			currentParentFitness = currentOffspringFitness;
		}

	L19:
		++lm;
		if (lm - currentConfiguration.n * currentConfiguration.lr < 0)
			goto L7;

		k = 1;
		if ((i_1 = le - l[0] - currentConfiguration.n - currentConfiguration.n) >= 0)
			if (i_1 == 0)
				--k;
			else
				k = k-2;
		i_1 = currentConfiguration.n;
		for (i = 1; i <= i_1; ++i)
		{
			// Computing MAX
			d_2 = currentConfiguration.sm[i] * pow(currentConfiguration.sn, k),
			d_3 = (d_1 = currentConfiguration.xb[i], doubleAbs(d_1)) * currentConfiguration.eb,
			d_2 = max(d_2,d_3);
			currentConfiguration.sm[i] = max(d_2, currentConfiguration.ea);
		}
		for (k = 1; k <= 9; ++k)
			l[k - 1] = l[k];
		l[9] = le;
		lm = 0L;
		++lc;
		if (lc - currentConfiguration.ls * 10L < 0L)
			goto L7;
		if ((fc - currentConfiguration.fb - currentConfiguration.ec > 0.) &&
		((fc - currentConfiguration.fb) / currentConfiguration.ed - doubleAbs(fc) > 0.))
		{
			lc = 0L;
			fc = currentConfiguration.fb;

			goto L7;
		}
		currentConfiguration.lf = iSign(c__2, currentConfiguration.lf);
		// Best set of values have been found
		return(COMPLETE);
}



// FUNCTION: doubleAbs
// AUTHOR: Andrew Volz
// CREATED: 24/08/96
// LAST MODIFIED: 24/08/96
// DESCRIPTION: Returns the absolute value of a double.

double EvolutionStrategy::doubleAbs(double number)
{
	return(number >= 0 ? number : -number);
}



// FUNCTION: iSign
// AUTHOR: Schwefel, modifications by Andrew Volz
// CREATED: Unknown
// LAST MODIFIED: Unknown

long int EvolutionStrategy::iSign(long int a, long int b)
{
	long int x;

	x = (a >= 0 ? a : -a);
	return(b >= 0 ? x : -x);
}



// FUNCTION: getParentFitness
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Returns the fitness of the current parent of the population.

double EvolutionStrategy::getParentFitness() const
{
	return(currentParentFitness);
}



// FUNCTION: getOffspringFitness
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Returns the fitness of the current offspring of the population.

double EvolutionStrategy::getOffspringFitness() const
{
	return(currentOffspringFitness);
}



// FUNCTION: copyConfiguration
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Copies the initial configuration to the current configuration.

int EvolutionStrategy::copyConfiguration(const ESConfiguration &initialConfiguration)
{
	currentConfiguration.n = initialConfiguration.n;
	currentConfiguration.m = initialConfiguration.m;
	currentConfiguration.lf = initialConfiguration.lf;
	currentConfiguration.lr = initialConfiguration.lr;
	currentConfiguration.ls = initialConfiguration.ls;
	currentConfiguration.tm = initialConfiguration.tm;
	currentConfiguration.ea = initialConfiguration.ea;
	currentConfiguration.eb = initialConfiguration.eb;
	currentConfiguration.ec = initialConfiguration.ec;
	currentConfiguration.ed = initialConfiguration.ed;
	currentConfiguration.sn = initialConfiguration.sn;
	currentConfiguration.fb = initialConfiguration.fb;
	copyParameters(currentConfiguration.xb, initialConfiguration.xb, initialConfiguration.n+1);
	copyParameters(currentConfiguration.sm, initialConfiguration.sm, initialConfiguration.n+1);
	return(TRUE);
}



// FUNCTION: uniformRandomNumber
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Returns a uniformly distributed random number.

double EvolutionStrategy::uniformRandomNumber()
{
	double number;
	number = (double) rand();
	number /= RAND_MAX;
	return(number);
}



// FUNCTION: getParentValues
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Gets the values of the current parent.

int EvolutionStrategy::getParentValues(double xVector[])
{
	copyParameters(xVector, currentConfiguration.xb, currentConfiguration.n+1);
	return(TRUE);
}



// FUNCTION: getOffspringValues
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Gets the values of the current offspring.

int EvolutionStrategy::getOffspringValues(double xVector[])
{
	copyParameters(xVector, currentConfiguration.x, currentConfiguration.n+1);
	return(TRUE);
}



// FUNCTION: getOffspringValues
// AUTHOR: Andrew Volz
// CREATED: 27/08/96
// LAST MODIFIED: 27/08/96
// DESCRIPTION: Gets a code representing the current status of the evolution strategy.

long int EvolutionStrategy::getStatusCode() const
{
	return(currentConfiguration.lf);
}



// FUNCTION: getNumberOfMutations
// AUTHOR: Andrew Volz
// CREATED: 29/08/96
// LAST MODIFIED: 29/08/96
// DESCRIPTION: Returns a value describing the number of mutations which have
//		currently been performed.

int EvolutionStrategy::getNumberOfMutations() const
{
	return(noOfMutations);
}



// FUNCTION: getNumberOfIterations
// AUTHOR: Andrew Volz
// CREATED: 07/09/96
// LAST MODIFIED: 07/09/96
// DESCRIPTION: Returns a value describing the number of iterations which have
//		currently been performed.

int EvolutionStrategy::getNumberOfIterations() const
{
	return(noOfIterations);
}