From 6f03c0b106d6b49b2d7aac026a19a735529f0573 Mon Sep 17 00:00:00 2001 From: swaits Date: Fri, 25 May 2007 01:24:26 +0000 Subject: [PATCH] use normal distribution on initial random individual creation initial default range tweaked way down rename "POP" template parameter to "POPSIZE" added RunOneGeneration() method - generations can now be run one at a time added SetRange(i,min,max) method added a decent amount of comments through main algorithm and class declaration --- de.h | 68 ++++++++------ de.inl | 210 ++++++++++++++++++++++++++------------------ derand.h | 2 + derand.inl | 10 +++ piecewiseapprox.cpp | 90 +++++++++++++++++++ 5 files changed, 267 insertions(+), 113 deletions(-) create mode 100644 piecewiseapprox.cpp diff --git a/de.h b/de.h index 93bdf8f..1d4c767 100644 --- a/de.h +++ b/de.h @@ -8,9 +8,8 @@ namespace DE { - - template + template class Engine { public: @@ -18,58 +17,69 @@ namespace DE Engine(); virtual ~Engine(); - void Reset(); - const double* GetBest() const; void SetRange(double minimum, double maximum); + void SetRange(unsigned int i, double minimum, double maximum); void SetRange(const Vector& minimum, const Vector& maximum); virtual double CalculateError(const double testcase[DIM], bool& stop) = 0; + void Reset(); bool Solve(unsigned int maxgenerations = 1000000000); + bool RunOneGeneration(); private: - void Select(unsigned int candidate, unsigned int *r1 = 0, unsigned int *r2 = 0, unsigned int *r3 = 0, unsigned int *r4 = 0, unsigned int *r5 = 0); + // randomly select unique individuals + void Select( + unsigned int candidate, + unsigned int *r1 = 0, + unsigned int *r2 = 0, + unsigned int *r3 = 0, + unsigned int *r4 = 0, + unsigned int *r5 = 0 + ); + // print a single individual to stdout void Dump(const double& fitness, Vector& individual) const; - void MakeTrial_best1exp(unsigned int candidate, Vector& trial); - void MakeTrial_rand1exp(unsigned int candidate, Vector& trial); - void MakeTrial_randtobest1exp(unsigned int candidate, Vector& trial); - void MakeTrial_best2exp(unsigned int candidate, Vector& trial); - void MakeTrial_rand2exp(unsigned int candidate, Vector& trial); - void MakeTrial_best1bin(unsigned int candidate, Vector& trial); - void MakeTrial_rand1bin(unsigned int candidate, Vector& trial); - void MakeTrial_randtobest1bin(unsigned int candidate, Vector& trial); - void MakeTrial_best2bin(unsigned int candidate, Vector& trial); - void MakeTrial_rand2bin(unsigned int candidate, Vector& trial); - - private: + // DE algorithm + void MakeTrial_best1exp (unsigned int candidate, Vector& trial); + void MakeTrial_rand1exp (unsigned int candidate, Vector& trial); + void MakeTrial_randtobest1exp (unsigned int candidate, Vector& trial); + void MakeTrial_best2exp (unsigned int candidate, Vector& trial); + void MakeTrial_rand2exp (unsigned int candidate, Vector& trial); + void MakeTrial_best1bin (unsigned int candidate, Vector& trial); + void MakeTrial_rand1bin (unsigned int candidate, Vector& trial); + void MakeTrial_randtobest1bin (unsigned int candidate, Vector& trial); + void MakeTrial_best2bin (unsigned int candidate, Vector& trial); + void MakeTrial_rand2bin (unsigned int candidate, Vector& trial); + // constants, parameters static const double DEFAULTSCALE; static const double DEFAULTCROSSOVER; static const double DEFAULTRANGE; static const double BIGDOUBLE; + double scale; + double crossover; + // entropy source + Random prng; - private: - - Random prng; - + // population Vector best; - Vector minimum; - Vector maximum; - Vector population[POP]; - double fitness[POP]; - - double scale; - double crossover; + Vector mean; + Vector variance; + Vector population[POPSIZE]; + double fitness[POPSIZE]; double bestfitness; - }; + // algorithm state + bool success; + unsigned int generation; + }; }; diff --git a/de.inl b/de.inl index e60735f..de8577d 100644 --- a/de.inl +++ b/de.inl @@ -10,43 +10,44 @@ #endif -template -const double DE::Engine::DEFAULTSCALE = 0.85; +template +const double DE::Engine::DEFAULTSCALE = 0.85; -template -const double DE::Engine::DEFAULTCROSSOVER = 1.0; +template +const double DE::Engine::DEFAULTCROSSOVER = 1.0; -template -const double DE::Engine::DEFAULTRANGE= 100.0; +template +const double DE::Engine::DEFAULTRANGE= 0.1; -template -const double DE::Engine::BIGDOUBLE = 1.79e308; // close to max double +template +const double DE::Engine::BIGDOUBLE = 1.79e308; // close to max double -template -inline DE::Engine::Engine() +template +inline DE::Engine::Engine() { + // set defaults scale = DEFAULTSCALE; crossover = DEFAULTCROSSOVER; - - for ( unsigned int i=0;iminimum[i] = -DEFAULTRANGE; - this->maximum[i] = -(this->minimum[i]); + this->mean[i] = 0.0; + this->variance[i] = DEFAULTRANGE; } + // create an initial population Reset(); return; } -template -inline DE::Engine::~Engine() +template +inline DE::Engine::~Engine() { return; } -template -inline void DE::Engine::Reset() +template +inline void DE::Engine::Reset() { // reset our best bestfitness = BIGDOUBLE; @@ -56,7 +57,7 @@ inline void DE::Engine::Reset() } // reset all vectors - for ( unsigned int i=0;i::Reset() // make a new individual for ( unsigned int j=0;j -inline void DE::Engine::SetRange(double minimum, double maximum) +template +inline void DE::Engine::SetRange(double minimum, double maximum) { for (unsigned int i=0;iminimum[i] = minimum; - this->maximum[i] = maximum; + SetRange(i,minimum,maximum); } return; } -template -inline void DE::Engine::SetRange(const Vector& minimum, const Vector& maximum) +template +inline void DE::Engine::SetRange(unsigned int i, double minimum, double maximum) { - this->minimum = minimum; - this->maximum = maximum; + if (imean[i] = (maximum+minimum)/2.0; + this->variance[i] = (maximum-minimum)/2.0; + } return; } -template -inline const double* DE::Engine::GetBest() const +template +inline void DE::Engine::SetRange(const Vector& minimum, const Vector& maximum) +{ + for (unsigned int i=0;i +inline const double* DE::Engine::GetBest() const { return best; } -template -inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1, unsigned int *r2, unsigned int *r3, unsigned int *r4, unsigned int *r5) +template +inline void DE::Engine::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,POP-1); + *r1 = prng.RandU32(0,POPSIZE-1); } while ( *r1 == candidate ); } @@ -114,7 +139,7 @@ inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1 { do { - *r2 = prng.RandU32(0,POP-1); + *r2 = prng.RandU32(0,POPSIZE-1); } while ( (*r2 == candidate) || (*r2 == *r1) ); } @@ -123,7 +148,7 @@ inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1 { do { - *r3 = prng.RandU32(0,POP-1); + *r3 = prng.RandU32(0,POPSIZE-1); } while ( (*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1) ); } @@ -132,7 +157,7 @@ inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1 { do { - *r4 = prng.RandU32(0,POP-1); + *r4 = prng.RandU32(0,POPSIZE-1); } while ( (*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1) ); } @@ -141,7 +166,7 @@ inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1 { do { - *r5 = prng.RandU32(0,POP-1); + *r5 = prng.RandU32(0,POPSIZE-1); } while ( (*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3) || (*r5 == *r2) || (*r5 == *r1) ); } @@ -149,49 +174,66 @@ inline void DE::Engine::Select(unsigned int candidate, unsigned int *r1 return; } -// TODO: change this to run one generation at a time??? -template -inline bool DE::Engine::Solve(unsigned int maxgenerations) +template +inline bool DE::Engine::Solve(unsigned int maxgenerations) { - bool success = false; - + // spawn a new generation Reset(); - for (unsigned int generation=0;generation +bool DE::Engine::RunOneGeneration() +{ + printf("====== GENERATION %07d =============================================\n\n",++generation); + + // attempt to improve each individual in the population + for (int candidate=0;candidate trial; + MakeTrial_randtobest1exp(candidate,trial); + + // determine fitness of trial individual + double trialfitness = CalculateError(trial,success); + + // see if the trial improved current candidate in population + if ( trialfitness < fitness[candidate] ) { - Vector trial; - MakeTrial_randtobest1exp(candidate,trial); + // it did improve, save it + fitness[candidate] = trialfitness; + population[candidate] = trial; - double trialfitness = CalculateError(trial,success); - - if ( trialfitness < fitness[candidate] ) + // see if it also improved upon our current best + if ( trialfitness < bestfitness ) { - fitness[candidate] = trialfitness; - population[candidate] = trial; + // it's better than the best, save it + bestfitness = trialfitness; + best = trial; - if ( trialfitness < bestfitness ) - { - bestfitness = trialfitness; - best = trial; - - Dump(bestfitness,best); - } + // output new best + Dump(bestfitness,best); } } } - printf("====== RUN COMPLETE =============================================\n\n"); - Dump(bestfitness,best); - return success; } -template -inline void DE::Engine::Dump(const double& fitness, DE::Vector& individual) const +template +inline void DE::Engine::Dump(const double& fitness, DE::Vector& individual) const { printf(" fitness = %24.16f\n",fitness); for (unsigned int i=0;i::Dump(const double& fitness, DE::Vector& in printf("\n"); } -template -inline void DE::Engine::MakeTrial_best1exp(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_best1exp(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2; unsigned int n; @@ -220,8 +262,8 @@ inline void DE::Engine::MakeTrial_best1exp(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_rand1exp(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_rand1exp(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3; unsigned int n; @@ -239,8 +281,8 @@ inline void DE::Engine::MakeTrial_rand1exp(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_randtobest1exp(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_randtobest1exp(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2; unsigned int n; @@ -258,8 +300,8 @@ inline void DE::Engine::MakeTrial_randtobest1exp(unsigned int candidate return; } -template -inline void DE::Engine::MakeTrial_best2exp(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_best2exp(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3, r4; unsigned int n; @@ -277,8 +319,8 @@ inline void DE::Engine::MakeTrial_best2exp(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_rand2exp(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_rand2exp(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3, r4, r5; unsigned int n; @@ -296,8 +338,8 @@ inline void DE::Engine::MakeTrial_rand2exp(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_best1bin(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_best1bin(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2; unsigned int n; @@ -318,8 +360,8 @@ inline void DE::Engine::MakeTrial_best1bin(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_rand1bin(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_rand1bin(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3; unsigned int n; @@ -340,8 +382,8 @@ inline void DE::Engine::MakeTrial_rand1bin(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_randtobest1bin(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_randtobest1bin(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2; unsigned int n; @@ -362,8 +404,8 @@ inline void DE::Engine::MakeTrial_randtobest1bin(unsigned int candidate return; } -template -inline void DE::Engine::MakeTrial_best2bin(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_best2bin(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3, r4; unsigned int n; @@ -384,8 +426,8 @@ inline void DE::Engine::MakeTrial_best2bin(unsigned int candidate, DE:: return; } -template -inline void DE::Engine::MakeTrial_rand2bin(unsigned int candidate, DE::Vector& trial) +template +inline void DE::Engine::MakeTrial_rand2bin(unsigned int candidate, DE::Vector& trial) { unsigned int r1, r2, r3, r4, r5; unsigned int n; diff --git a/derand.h b/derand.h index e312bfa..04db5e2 100644 --- a/derand.h +++ b/derand.h @@ -25,6 +25,8 @@ namespace DE double RandDouble(); double RandDouble(double min, double max); + + double RandGaussian(double mean, double variance); u32 RandU32(); u32 RandU32(u32 min, u32 max); diff --git a/derand.inl b/derand.inl index e6db650..7eac604 100644 --- a/derand.inl +++ b/derand.inl @@ -50,6 +50,16 @@ inline double DE::Random::RandDouble(double min, double max) return ( RandDouble() * (max - min) + min ); } +inline double DE::Random::RandGaussian(double mean, double variance) +{ + double sum = 0.0; + for (int i=0;i<12;++i) + { + sum += RandDouble(); + } + return mean + ( variance * ( sum - 6.0 ) ); +} + inline DE::u32 DE::Random::RandU32() { return Rand(); diff --git a/piecewiseapprox.cpp b/piecewiseapprox.cpp new file mode 100644 index 0000000..f4b77ee --- /dev/null +++ b/piecewiseapprox.cpp @@ -0,0 +1,90 @@ +#include +#include + +#include "de.h" + +// the order of the polynomial (actually the number of coefficients, i.e. 1 = constant, 2 = line, etc.) +#define NUM_COEFFICIENTS 7 + +// the function we're searching for +double TargetFunction(double x) +{ + //return 1.0 / (1.0 + exp(-4.9*x)); // Ken Stanley's NEAT activation + return (1.0 / (1.0 + exp(-4.9*x)) - 0.5) * 2.0; // " " (signed version) + //return (1.0 / (1.0 + exp(-x)) - 0.5) * 2.0; // signed sigmoid + //return exp(-x*x); // gaussian + //return tanh(x); // tanh + //return (2.0/(1.0+exp(-x))-1.0); // tanh approximation + //return (x<0.0 ? -1.0 : 1.0) * sqrt(abs(x)); // signed root + //return sin(x); // sin +} + + +// a helper function used in calculating the polynomials +double CalculatePiecewise(double x, int n_coeffs, const double coeffs[]) +{ + double abs_x = (x < 0.0 ? -x : x); + const double* pcoeffs; + if ( abs_x < coeffs[6] ) + pcoeffs = &coeffs[0]; + else + pcoeffs = &coeffs[3]; + + return (pcoeffs[0] + pcoeffs[1]*abs_x + pcoeffs[2]*abs_x*abs_x) * (x < 0.0 ? -1.0 : 1.0); +} + + +// our DE implementation +class PolySearch: public DE::Engine +{ + public: + // this method is called by the DE::Engine::Solve() to test a single DE individual + double CalculateError(const double testcase[NUM_COEFFICIENTS], bool& stop) + { + double error = 0.0; + unsigned int test_count = 0; + + // accumulate squared error + for ( double i=-2.5;i<=2.5;i+=0.01 ) + { + // test + double test = CalculatePiecewise(i,NUM_COEFFICIENTS,testcase); + double actual = TargetFunction(i); + + // accumulate error^2 + double e = test - actual; + error += e * e; + + // inc test count (divisor to get from SE to MSE) + test_count++; + } + + // convert SE to RMSE + error = sqrt(error / (double)test_count); + + // set stop flag if the error is really low + if ( error < 0.0000001 ) + { + stop = true; + } + + // done + return error; + } +}; + + +// main program +int main() +{ + PolySearch de; + + // setup, adjust from default range + //de.SetRange(-0.1,0.1); + + // reset, and run the DE + de.Solve(1000); + + // done + return 0; +}