/*
The MIT License

Copyright (c) 2010, Stephen Waits <steve@waits.net>

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/



#ifndef __de_inl__
#define __de_inl__


#include <cstdio>
#include <cmath>


#if !(defined(__de_h__))
#	 error "de.inl should only be included by de.h"
#endif


template <unsigned int DIM, unsigned int POPSIZE>
const double DE::Engine<DIM,POPSIZE>::DEFAULTSCALE = 0.85;

template <unsigned int DIM, unsigned int POPSIZE>
const double DE::Engine<DIM,POPSIZE>::DEFAULTCROSSOVER = 1.0;

template <unsigned int DIM, unsigned int POPSIZE>
const double DE::Engine<DIM,POPSIZE>::DEFAULTRANGE= 0.1;

template <unsigned int DIM, unsigned int POPSIZE>
const double DE::Engine<DIM,POPSIZE>::BIGDOUBLE = 1.79e308; // close to max double

template <unsigned int DIM, unsigned int POPSIZE>
inline DE::Engine<DIM,POPSIZE>::Engine()
{
	// set defaults
	scale     = DEFAULTSCALE;
	crossover = DEFAULTCROSSOVER;
	for ( unsigned int i=0;i<DIM;++i )
	{
		this->mean[i]     = 0.0;
		this->variance[i] = DEFAULTRANGE;
	}

	// create an initial population
	Reset();

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline DE::Engine<DIM,POPSIZE>::~Engine()
{
	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::Reset()
{
	// reset our best
	bestfitness = BIGDOUBLE;
	for ( unsigned int i=0;i<DIM;i++ )
	{
		best[i] = 0.0;
	}

	// reset all vectors
	for ( unsigned int i=0;i<POPSIZE;i++ )
	{
		// energy tied for worst
		fitness[i] = bestfitness;

		// make a new individual
		for ( unsigned int j=0;j<DIM;j++ )
		{
			population[i][j] = prng.RandGaussian(mean[j],variance[j]);
		}
	}

	// reset state
	generation = 0;
	success    = false;

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::SetRange(double minimum, double maximum)
{
	for (unsigned int i=0;i<DIM;i++)
	{
		SetRange(i,minimum,maximum);
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::SetRange(unsigned int i, double minimum, double maximum)
{
	if (i<DIM)
	{
		this->mean[i]     = (maximum+minimum)/2.0;
		this->variance[i] = (maximum-minimum)/2.0;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::SetRange(const Vector<DIM>& minimum, const Vector<DIM>& maximum)
{
	for (unsigned int i=0;i<DIM;i++)
	{
		SetRange(i,minimum[i],maximum[i]);
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline const double* DE::Engine<DIM,POPSIZE>::GetBest() const
{
	return best;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::Select(
	unsigned int candidate,
	unsigned int *r1,
	unsigned int *r2,
	unsigned int *r3,
	unsigned int *r4,
	unsigned int *r5
)
{
	if ( r1 )
	{
		do
		{
			*r1 = prng.RandU32(0,POPSIZE-1);
		}
		while ( *r1 == candidate );
	}

	if ( r2 )
	{
		do
		{
			*r2 = prng.RandU32(0,POPSIZE-1);
		}
		while ( (*r2 == candidate) || (*r2 == *r1) );
	}

	if ( r3 )
	{
		do
		{
			*r3 = prng.RandU32(0,POPSIZE-1);
		}
		while ( (*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1) );
	}

	if ( r4 )
	{
		do
		{
			*r4 = prng.RandU32(0,POPSIZE-1);
		}
		while ( (*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1) );
	}

	if ( r5 )
	{
		do
		{
			*r5 = prng.RandU32(0,POPSIZE-1);
		}
		while ( (*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3) || (*r5 == *r2) || (*r5 == *r1) );
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline bool DE::Engine<DIM,POPSIZE>::Solve(unsigned int maxgenerations)
{
	// spawn a new generation
	Reset();

	// run generations until max or success
	while(!success && generation<maxgenerations)
	{
		success = RunOneGeneration();
	}

	// output
	printf("====== RUN COMPLETE       =============================================\n\n");
	Dump(bestfitness,best);

	// finished
	return success;
}

template <unsigned int DIM, unsigned int POPSIZE>
bool DE::Engine<DIM,POPSIZE>::RunOneGeneration()
{
	printf("====== GENERATION %07d =============================================\n\n",++generation);

	// attempt to improve each individual in the population
	for (int candidate=0;candidate<POPSIZE && !success;candidate++)
	{
		// create a new trial individual to test
		Vector<DIM> trial;
		MakeTrial_randtobest1exp(candidate,trial);

		// determine fitness of trial individual
		double trialfitness = TestFitness(trial,success);

		// see if the trial improved current candidate in population
		if ( trialfitness < fitness[candidate] )
		{
			// it did improve, save it
			fitness[candidate]    = trialfitness;
			population[candidate] = trial;

			// see if it also improved upon our current best
			if ( trialfitness < bestfitness )
			{
				// it's better than the best, save it
				bestfitness = trialfitness;
				best        = trial;

				// output new best
				Dump(bestfitness,best);
			}
		}
	}

	return success;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::Dump(const double& fitness, DE::Vector<DIM>& individual) const
{
	printf(" fitness = %24.16f\n",fitness);
	for (unsigned int i=0;i<DIM;i++)
	{
		printf("   [%03d] = %24.16e\n",i,individual[i]);
	}
	printf("\n");
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_best1exp(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2;
	unsigned int n;

	Select(candidate, &r1, &r2);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; (prng.RandDouble(0.0,1.0) < crossover) && (i < DIM); i++) 
	{
		trial[n] = best[n] + scale * (population[r1][n] - population[r2][n]);
		n        = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_rand1exp(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; (prng.RandDouble(0.0,1.0) <= crossover) && (i < DIM); i++) 
	{
		trial[n] = population[r1][n] + scale * (population[r2][n] - population[r3][n]);
		n        = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_randtobest1exp(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2;
	unsigned int n;

	Select(candidate, &r1, &r2);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; (prng.RandDouble(0.0,1.0) <= crossover) && (i < DIM); i++) 
	{
		trial[n] += scale * (best[n] - trial[n]) + scale * (population[r1][n] - population[r2][n]);
		n         = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_best2exp(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3, r4;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3, &r4);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; (prng.RandDouble(0.0,1.0) <= crossover) && (i < DIM); i++) 
	{
		trial[n] = best[n] + scale * (population[r1][n] + population[r2][n] - population[r3][n] - population[r4][n]);
		n        = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_rand2exp(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3, r4, r5;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3 ,&r4, &r5);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; (prng.RandDouble(0.0,1.0) <= crossover) && (i < DIM); i++) 
	{
		trial[n] = population[r1][n] + scale * (population[r2][n] + population[r3][n] - population[r4][n] - population[r5][n]);
		n        = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_best1bin(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2;
	unsigned int n;

	Select(candidate, &r1, &r2);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; i < DIM; i++) 
	{
		if ( (prng.RandDouble(0.0,1.0) <= crossover) || (i == (DIM - 1)) )
		{
			trial[n] = best[n] + scale * (population[r1][n] - population[r2][n]);
		}
		n = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_rand1bin(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; i < DIM; i++) 
	{
		if ( (prng.RandDouble(0.0,1.0) <= crossover) || (i == (DIM - 1)) )
		{
			trial[n] = population[r1][n] + scale * (population[r2][n] - population[r3][n]);
		}
		n = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_randtobest1bin(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2;
	unsigned int n;

	Select(candidate, &r1, &r2);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; i < DIM; i++) 
	{
		if ( (prng.RandDouble(0.0,1.0) <= crossover) || (i == (DIM - 1)) )
		{
			trial[n] += scale * (best[n] - trial[n]) + scale * (population[r1][n] - population[r2][n]);
		}
		n = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_best2bin(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3, r4;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3, &r4);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; i < DIM; i++) 
	{
		if ( (prng.RandDouble(0.0,1.0) <= crossover) || (i == (DIM - 1)) )
		{
			trial[n] = best[n] + scale * (population[r1][n] + population[r2][n] - population[r3][n] - population[r4][n]);
		}
		n = (n + 1) % DIM;
	}

	return;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::MakeTrial_rand2bin(unsigned int candidate, DE::Vector<DIM>& trial)
{
	unsigned int r1, r2, r3, r4, r5;
	unsigned int n;

	Select(candidate, &r1, &r2, &r3, &r4, &r5);
	n = prng.RandU32(0,DIM-1);

	trial = population[candidate];
	for (unsigned int i=0; i < DIM; i++) 
	{
		if ( (prng.RandDouble(0.0,1.0) <= crossover) || (i == (DIM - 1)) )
		{
			trial[n] = population[r1][n] + scale * (population[r2][n] + population[r3][n] - population[r4][n] - population[r5][n]);
		}
		n = (n + 1) % DIM;
	}

	return;
}




template <unsigned int DIM, unsigned int POPSIZE>
inline DE::Engine<DIM,POPSIZE>::ErrorAccumulator::ErrorAccumulator()
{
	Clear();
}

template <unsigned int DIM, unsigned int POPSIZE>
inline DE::Engine<DIM,POPSIZE>::ErrorAccumulator::~ErrorAccumulator()
{
	// empty
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::ErrorAccumulator::Clear()
{
	num_tests         = 0;
	sum_error_squared = 0.0;
	error_rmse        = 0.0; // initial error == zero
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::ErrorAccumulator::AddTestCase(double expected, double actual)
{
	double e           = expected - actual; // compute error
	sum_error_squared += e * e;             // add to sum of error^2
	num_tests         += 1;                 // increment total test count
	error_rmse         = -1.0;              // mark as dirty
}

template <unsigned int DIM, unsigned int POPSIZE>
inline double DE::Engine<DIM,POPSIZE>::ErrorAccumulator::GetRMSE()
{
	if ( error_rmse < 0.0 ) // recalculate if set to "dirty"
		CalculateError();

	return error_rmse;
}

template <unsigned int DIM, unsigned int POPSIZE>
inline void DE::Engine<DIM,POPSIZE>::ErrorAccumulator::CalculateError()
{
	// convert squared error to root mean squared error
	error_rmse = ::sqrt(sum_error_squared / (double)num_tests);
}



#endif // __de_inl__