From 70650ba7531b83e3bb899cbcbc99cff4e319c47d Mon Sep 17 00:00:00 2001 From: swaits Date: Fri, 20 May 2005 18:10:11 +0000 Subject: [PATCH] initial import of differential evolution engine --- de.h | 80 +++++++++++ de.inl | 407 +++++++++++++++++++++++++++++++++++++++++++++++++++++ derand.cpp | 1 + derand.h | 76 ++++++++++ derand.inl | 69 +++++++++ devec.h | 63 +++++++++ devec.inl | 142 +++++++++++++++++++ 7 files changed, 838 insertions(+) create mode 100644 de.h create mode 100644 de.inl create mode 100644 derand.cpp create mode 100644 derand.h create mode 100644 derand.inl create mode 100644 devec.h create mode 100644 devec.inl diff --git a/de.h b/de.h new file mode 100644 index 0000000..1f2648a --- /dev/null +++ b/de.h @@ -0,0 +1,80 @@ +#ifndef __de_h__ +#define __de_h__ + + +#include "devec.h" +#include "derand.h" + + +namespace DE +{ + + + template + class Engine + { + public: + + Engine(); + virtual ~Engine(); + + void Reset(); + + const double* GetBest() const; + + void SetRange(double minimum, double maximum); + void SetRange(const Vector& minimum, const Vector& maximum); + + virtual double CalculateError(const double testcase[DIM], bool& stop) = 0; + + bool Solve(unsigned int maxgenerations = 1000000000); + + 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); + + 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: + + static const double DEFAULTSCALE; + static const double DEFAULTCROSSOVER; + static const double BIGDOUBLE; + + + private: + + Random prng; + + Vector best; + Vector minimum; + Vector maximum; + Vector population[POP]; + double fitness[POP]; + + double scale; + double crossover; + double bestfitness; + + }; + + +}; + + +#include "de.inl" + + +#endif // __de_h__ + diff --git a/de.inl b/de.inl new file mode 100644 index 0000000..9ec5b81 --- /dev/null +++ b/de.inl @@ -0,0 +1,407 @@ +#ifndef __de_inl__ +#define __de_inl__ + + +#include + + +#if !(defined(__de_h__)) +# error "de.inl should only be included by de.h" +#endif + + +template +const double DE::Engine::DEFAULTSCALE = 0.85; + +template +const double DE::Engine::DEFAULTCROSSOVER = 1.0; + +template +const double DE::Engine::BIGDOUBLE = 1.0e30; + +template +inline DE::Engine::Engine() +{ + scale = DEFAULTSCALE; + crossover = DEFAULTCROSSOVER; + + for ( unsigned int i=0;iminimum[i] = -BIGDOUBLE; + this->maximum[i] = -(this->minimum[i]); + } + + Reset(); + + return; +} + +template +inline DE::Engine::~Engine() +{ + return; +} + +template +inline void DE::Engine::Reset() +{ + // reset our best + bestfitness = BIGDOUBLE; + for ( unsigned int i=0;i +inline void DE::Engine::SetRange(double minimum, double maximum) +{ + for (unsigned int i=0;iminimum[i] = minimum; + this->maximum[i] = maximum; + } + + return; +} + +template +inline void DE::Engine::SetRange(const Vector& minimum, const Vector& maximum) +{ + this->minimum = minimum; + this->maximum = maximum; + + return; +} + +template +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) +{ + if ( r1 ) + { + do + { + *r1 = prng.RandU32(0,POP-1); + } + while ( *r1 == candidate ); + } + + if ( r2 ) + { + do + { + *r2 = prng.RandU32(0,POP-1); + } + while ( (*r2 == candidate) || (*r2 == *r1) ); + } + + if ( r3 ) + { + do + { + *r3 = prng.RandU32(0,POP-1); + } + while ( (*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1) ); + } + + if ( r4 ) + { + do + { + *r4 = prng.RandU32(0,POP-1); + } + while ( (*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1) ); + } + + if ( r5 ) + { + do + { + *r5 = prng.RandU32(0,POP-1); + } + while ( (*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3) || (*r5 == *r2) || (*r5 == *r1) ); + } + + return; +} + +// TODO: change this to run one generation at a time??? +template +inline bool DE::Engine::Solve(unsigned int maxgenerations) +{ + bool success = false; + + Reset(); + + for (unsigned int generation=0;generation trial; + MakeTrial_randtobest1bin(candidate,trial); + + double trialfitness = CalculateError(trial,success); + + if ( trialfitness < fitness[candidate] ) + { + fitness[candidate] = trialfitness; + population[candidate] = trial; + + if ( trialfitness < bestfitness ) + { + bestfitness = trialfitness; + best = trial; + + 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 +{ + printf(" fitness = %60.30f\n",fitness); + for (unsigned int i=0;i +inline void DE::Engine::MakeTrial_best1exp(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_rand1exp(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_randtobest1exp(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_best2exp(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_rand2exp(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_best1bin(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_rand1bin(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_randtobest1bin(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_best2bin(unsigned int candidate, DE::Vector& 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 +inline void DE::Engine::MakeTrial_rand2bin(unsigned int candidate, DE::Vector& 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; +} + +#endif // __de_inl__ + diff --git a/derand.cpp b/derand.cpp new file mode 100644 index 0000000..d3f5a12 --- /dev/null +++ b/derand.cpp @@ -0,0 +1 @@ + diff --git a/derand.h b/derand.h new file mode 100644 index 0000000..11a7fe2 --- /dev/null +++ b/derand.h @@ -0,0 +1,76 @@ +#ifndef __derand_h__ +#define __derand_h__ + +namespace DE +{ + + + typedef unsigned int u32; + typedef signed int s32; + + + class Random + { + public: + + void Seed(u32 seed); + + u32 Rand(); + u32 operator()(); + + float RandFloat(); + float RandFloat(float min, float max); + + double RandDouble(); + double RandDouble(double min, double max); + + u32 RandU32(); + u32 RandU32(u32 min, u32 max); + + s32 RandS32(); + s32 RandS32(s32 min, s32 max); + + + + private: + + static const u32 U32_MAX = (u32)0xffffffff; + + /** + * LCG multiplier by M. Lavaux and F. Janssens + */ + static const u32 MULTIPLIER = 1664525UL; + + /** + * LCG increment by M. Lavaux and F. Janssens + */ + static const u32 INCREMENT = 1013904223UL; + + /** + * LCG modulus set to 2^32, but not used due to MODULUS_ASSUMED + * above (since we're using a 32 bit word here) + * + * Set this to 0 if the modulus falls out automatically + * due to your data type (i.e. "addition modulo". For + * instance, if you're using a 32 bit integer, and your + * modulus is 2^32, then there's no need to perform the + * modulo since a == a % (2^32). (so the modulo is "free") + * + */ + static const u32 MODULUS = 0xffffffff; + + /** + * LCG state + */ + u32 seed; + }; + + +}; + + +#include "derand.inl" + + +#endif // __derand_h__ + diff --git a/derand.inl b/derand.inl new file mode 100644 index 0000000..764ee62 --- /dev/null +++ b/derand.inl @@ -0,0 +1,69 @@ +#ifndef __derand_inl__ +#define __derand_inl__ + + +#if !(defined(__derand_h__)) +# error "derand.inl should only be included by derand.h" +#endif + + +inline void DE::Random::Seed(u32 seed) +{ + this->seed = seed; +} + +inline DE::u32 DE::Random::Rand() +{ + seed = ( (seed * DE::Random::MULTIPLIER) + DE::Random::INCREMENT ); + seed %= DE::Random::MODULUS; + + return seed; +} + +inline DE::u32 DE::Random::operator()() +{ + return this->Rand(); +} + +inline float DE::Random::RandFloat() +{ + return (float)Rand() / (float)DE::Random::U32_MAX; +} + +inline float DE::Random::RandFloat(float min, float max) +{ + return ( RandFloat() * (max - min) + min ); +} + +inline double DE::Random::RandDouble() +{ + return (double)Rand() / (double)DE::Random::U32_MAX; +} + +inline double DE::Random::RandDouble(double min, double max) +{ + return ( RandDouble() * (max - min) + min ); +} + +inline DE::u32 DE::Random::RandU32() +{ + return Rand(); +} + +inline DE::u32 DE::Random::RandU32(DE::u32 min, DE::u32 max) +{ + return ( ( Rand() % (max - min) ) + min ); +} + +inline DE::s32 DE::Random::RandS32() +{ + return (DE::s32)Rand(); +} + +inline DE::s32 DE::Random::RandS32(DE::s32 min, DE::s32 max) +{ + return (DE::s32)( ( Rand() % (max - min) ) + min ); +} + +#endif // __derand_inl__ + diff --git a/devec.h b/devec.h new file mode 100644 index 0000000..e80f064 --- /dev/null +++ b/devec.h @@ -0,0 +1,63 @@ +#ifndef __devec_h__ +#define __devec_h__ + +namespace DE +{ + + + template + class Vector + { + public: + + // default constructor + Vector(); + + // destructor + ~Vector(); + + // copy constructor + Vector(const Vector& other); + + // assignment + Vector& operator= (const Vector& other); + + // addition + Vector& operator+= (const Vector& other); + Vector operator+ (const Vector& other); + + // subtraction + Vector& operator-= (const Vector& other); + Vector operator- (const Vector& other); + + // multiply by scalar + Vector& operator*= (const double& multiplier); + Vector operator* (const double& multiplier); + + // divide by scalar + Vector& operator/= (const double& divisor); + Vector operator/ (const double& divisor); + + // element access + const double& operator[] (unsigned int i) const; + double& operator[] (unsigned int i); + + // data access + operator const double* () const; + operator double* (); + + private: + + // data lives here + double tuple[DIM]; + }; + + +}; + + +#include "devec.inl" + + +#endif // __devec_h__ + diff --git a/devec.inl b/devec.inl new file mode 100644 index 0000000..c7c59b4 --- /dev/null +++ b/devec.inl @@ -0,0 +1,142 @@ +#ifndef __devec_inl__ +#define __devec_inl__ + + +#if !(defined(__devec_h__)) +# error "devec.inl should only be included by devec.h" +#endif + + +template +inline DE::Vector::Vector() +{ + for (unsigned int i=0;i +inline DE::Vector::~Vector() +{ + // empty + return; +} + +template +inline DE::Vector::Vector(const Vector& other) +{ + *this = other; + return; +} + +template +inline DE::Vector& DE::Vector::operator= (const Vector& other) +{ + for (unsigned int i=0;i +inline DE::Vector& DE::Vector::operator+= (const Vector& other) +{ + for (unsigned int i=0;i +inline DE::Vector DE::Vector::operator+ (const Vector& other) +{ + Vector tmp(*this); + tmp += other; + return tmp; +} + +template +inline DE::Vector& DE::Vector::operator-= (const Vector& other) +{ + for (unsigned int i=0;i +inline DE::Vector DE::Vector::operator- (const Vector& other) +{ + Vector tmp(*this); + tmp -= other; + return tmp; +} + +template +inline DE::Vector& DE::Vector::operator*= (const double& multiplier) +{ + for (unsigned int i=0;i +inline DE::Vector DE::Vector::operator* (const double& multiplier) +{ + Vector tmp(*this); + tmp *= multiplier; + return tmp; +} + +template +inline DE::Vector& DE::Vector::operator/= (const double& divisor) +{ + for (unsigned int i=0;i +inline DE::Vector DE::Vector::operator/ (const double& divisor) +{ + Vector tmp(*this); + tmp /= divisor; + return tmp; +} + +template +inline const double& DE::Vector::operator[] (unsigned int i) const +{ + return tuple[i]; +} + +template +inline double& DE::Vector::operator[] (unsigned int i) +{ + return tuple[i]; +} + +template +inline DE::Vector::operator const double* () const +{ + return tuple; +} + +template +inline DE::Vector::operator double* () +{ + return tuple; +} + + +#endif // __devec_inl__ +