add github mirror

add README.md with description from previous Bitbucket repo
add quaternion search example
2024-01-29 10:42:37 -07:00 · 2020-04-24 08:20:32 -07:00 · 2011-01-24 14:18:30 -08:00 · 2010-09-30 16:39:40 -07:00 · 2010-09-30 16:19:03 -07:00 · 2010-09-30 15:53:16 -07:00
14 changed files with 1266 additions and 130 deletions
@@ -0,0 +1,11 @@
+image: alpine/edge
+secrets:
+  - 88cce185-d2ec-4573-9546-e33f36f79bbf
+sources:
+  - git+ssh://git@git.sr.ht/~swaits/differential-evolution
+tasks:
+  - mirror-to-github: |
+      cd ~/differential-evolution
+      ssh-keyscan -t rsa github.com >> ~/.ssh/known_hosts
+      git remote add github git@github.com:swaits/differential-evolution.git
+      git push --mirror github
@@ -0,0 +1,21 @@
+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.
@@ -0,0 +1,2 @@
+Differential Evolution real stochastic optimization library, in C++, with
+examples.
@@ -1,3 +1,48 @@
+/*
+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.
+*/
+
+
+/*
+ * This is an implementation of Differential Evolution in C++. It's primarily
+ * inline and has no external dependencies.
+ *
+ * DE is a simple, but useful stochastic optimization algorithm created by
+ * Kenneth Price and Rainer Storn. This code is based, in part, on their 
+ * original C version.
+ *
+ * See http://www.icsi.berkeley.edu/~storn/code.html for more information.
+ *
+ * The minimum you need to use this is to:
+ * 1. inherit DE::Engine
+ * 2. implement a concrete version of TestFitness in your derived class
+ * 3. instance your class and run Solve
+ *
+ * See the examples/ directory for some semi-practical examples on how you
+ * might make use of this code.
+ *
+ */
+
 #ifndef __de_h__
 #define __de_h__

@@ -8,9 +53,8 @@

 namespace DE
 {
-
 	
-	template <unsigned int DIM, unsigned int POP=(DIM*10)>
+	template <unsigned int DIM, unsigned int POPSIZE=(DIM*10)>
 	class Engine
 	{
 	public:
@@ -18,58 +62,92 @@ 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<DIM>& minimum, const Vector<DIM>& maximum);
 		
-		virtual double CalculateError(const double testcase[DIM], bool& stop) = 0;
+		virtual double TestFitness(const double testcase[DIM], bool& stop) = 0;

+		void Reset();
 		bool Solve(unsigned int maxgenerations = 1000000000);
+		bool RunOneGeneration();
+
+	public:
+
+		// helper class for accumulating error and calculating
+		// root mean square error
+		class ErrorAccumulator
+		{
+		public:
+			ErrorAccumulator();
+			virtual ~ErrorAccumulator();
+
+			void Clear();
+			void AddTestCase(double expected, double actual);
+			double GetRMSE();
+
+		private:
+
+			void CalculateError();
+
+			unsigned int num_tests;         // # of test cases
+			double       sum_error_squared; // sum of error^2
+			double       error_rmse;        // root mean square error (-1 == "dirty" flag)
+		};

 	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<DIM>& individual) const;

-		void MakeTrial_best1exp(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_rand1exp(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_randtobest1exp(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_best2exp(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_rand2exp(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_best1bin(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_rand1bin(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_randtobest1bin(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_best2bin(unsigned int candidate, Vector<DIM>& trial);
-		void MakeTrial_rand2bin(unsigned int candidate, Vector<DIM>& trial);
-
-	private:
+		// DE algorithm
+		void MakeTrial_best1exp       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_rand1exp       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_randtobest1exp (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_best2exp       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_rand2exp       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_best1bin       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_rand1bin       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_randtobest1bin (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_best2bin       (unsigned int candidate, Vector<DIM>& trial);
+		void MakeTrial_rand2bin       (unsigned int candidate, Vector<DIM>& 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<DIM> best;
-		Vector<DIM> minimum;
-		Vector<DIM> maximum;
-		Vector<DIM> population[POP];
-		double      fitness[POP];
-
-		double      scale;
-		double      crossover;
+		Vector<DIM> mean;
+		Vector<DIM> variance;
+		Vector<DIM> population[POPSIZE];
+		double      fitness[POPSIZE];
 		double      bestfitness;

-	};
+		// algorithm state
+		bool         success;
+		unsigned int generation;

+	};

 };

@@ -1,8 +1,35 @@
+/*
+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__))
@@ -10,43 +37,44 @@
 #endif


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

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

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

-template <unsigned int DIM, unsigned int POP>
-const double DE::Engine<DIM,POP>::BIGDOUBLE = 1.79e308; // close to max double
+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 POP>
-inline DE::Engine<DIM,POP>::Engine()
+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++ )
+	for ( unsigned int i=0;i<DIM;++i )
 	{
-		this->minimum[i] = -DEFAULTRANGE;
-		this->maximum[i] = -(this->minimum[i]);
+		this->mean[i]     = 0.0;
+		this->variance[i] = DEFAULTRANGE;
 	}

+	// create an initial population
 	Reset();

 	return;
 }

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

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::Reset()
+template <unsigned int DIM, unsigned int POPSIZE>
+inline void DE::Engine<DIM,POPSIZE>::Reset()
 {
 	// reset our best
 	bestfitness = BIGDOUBLE;
@@ -56,7 +84,7 @@ inline void DE::Engine<DIM,POP>::Reset()
 	}

 	// reset all vectors
-	for ( unsigned int i=0;i<POP;i++ )
+	for ( unsigned int i=0;i<POPSIZE;i++ )
 	{
 		// energy tied for worst
 		fitness[i] = bestfitness;
@@ -64,48 +92,72 @@ inline void DE::Engine<DIM,POP>::Reset()
 		// make a new individual
 		for ( unsigned int j=0;j<DIM;j++ )
 		{
-			population[i][j] = prng.RandDouble(minimum[j],maximum[j]);
+			population[i][j] = prng.RandGaussian(mean[j],variance[j]);
 		}
 	}

+	// reset state
+	generation = 0;
+	success    = false;
+
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::SetRange(double minimum, double maximum)
+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++)
 	{
-		this->minimum[i] = minimum;
-		this->maximum[i] = maximum;
+		SetRange(i,minimum,maximum);
 	}

 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::SetRange(const Vector<DIM>& minimum, const Vector<DIM>& maximum)
+template <unsigned int DIM, unsigned int POPSIZE>
+inline void DE::Engine<DIM,POPSIZE>::SetRange(unsigned int i, double minimum, double maximum)
 {
-	this->minimum = minimum;
-	this->maximum = 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 POP>
-inline const double* DE::Engine<DIM,POP>::GetBest() const
+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 POP>
-inline void DE::Engine<DIM,POP>::Select(unsigned int candidate, unsigned int *r1, unsigned int *r2, unsigned int *r3, unsigned int *r4, unsigned int *r5)
+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,POP-1);
+			*r1 = prng.RandU32(0,POPSIZE-1);
 		}
 		while ( *r1 == candidate );
 	}
@@ -114,7 +166,7 @@ inline void DE::Engine<DIM,POP>::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 +175,7 @@ inline void DE::Engine<DIM,POP>::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 +184,7 @@ inline void DE::Engine<DIM,POP>::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 +193,7 @@ inline void DE::Engine<DIM,POP>::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 +201,66 @@ inline void DE::Engine<DIM,POP>::Select(unsigned int candidate, unsigned int *r1
 	return;
 }

-// TODO: change this to run one generation at a time???
-template <unsigned int DIM, unsigned int POP>
-inline bool DE::Engine<DIM,POP>::Solve(unsigned int maxgenerations)
+template <unsigned int DIM, unsigned int POPSIZE>
+inline bool DE::Engine<DIM,POPSIZE>::Solve(unsigned int maxgenerations)
 {
-	bool success = false;
-
+	// spawn a new generation
 	Reset();

-	for (unsigned int generation=0;generation<maxgenerations && !success;generation++)
+	// run generations until max or success
+	while(!success && generation<maxgenerations)
 	{
-		printf("====== GENERATION %07d =============================================\n\n",generation);
+		success = RunOneGeneration();
+	}

-		for (int candidate=0;candidate<POP && !success;candidate++)
+	// 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] )
 		{
-			Vector<DIM> 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 <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::Dump(const double& fitness, DE::Vector<DIM>& individual) const
+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++)
@@ -201,8 +270,8 @@ inline void DE::Engine<DIM,POP>::Dump(const double& fitness, DE::Vector<DIM>& in
 	printf("\n");
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_best1exp(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -220,8 +289,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_best1exp(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_rand1exp(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -239,8 +308,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_rand1exp(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_randtobest1exp(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -258,8 +327,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_randtobest1exp(unsigned int candidate
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_best2exp(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -277,8 +346,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_best2exp(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_rand2exp(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -296,8 +365,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_rand2exp(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_best1bin(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -318,8 +387,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_best1bin(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_rand1bin(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -340,8 +409,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_rand1bin(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_randtobest1bin(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -362,8 +431,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_randtobest1bin(unsigned int candidate
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_best2bin(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -384,8 +453,8 @@ inline void DE::Engine<DIM,POP>::MakeTrial_best2bin(unsigned int candidate, DE::
 	return;
 }

-template <unsigned int DIM, unsigned int POP>
-inline void DE::Engine<DIM,POP>::MakeTrial_rand2bin(unsigned int candidate, DE::Vector<DIM>& trial)
+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;
@@ -406,5 +475,55 @@ inline void DE::Engine<DIM,POP>::MakeTrial_rand2bin(unsigned int candidate, DE::
 	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__

@@ -1,3 +1,29 @@
+/*
+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 __derand_h__
 #define __derand_h__

@@ -25,6 +51,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);
@@ -1,3 +1,29 @@
+/*
+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 __derand_inl__
 #define __derand_inl__

@@ -50,6 +76,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();
@@ -1,3 +1,29 @@
+/*
+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 __devec_h__
 #define __devec_h__

@@ -1,3 +1,29 @@
+/*
+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 __devec_inl__
 #define __devec_inl__

@@ -0,0 +1,82 @@
+#include <cstdio>
+#include <cmath>
+
+#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<NUM_COEFFICIENTS>
+{
+	public:
+		// this method is called by the DE::Engine::Solve() to test a single DE individual
+		double TestFitness(const double testcase[NUM_COEFFICIENTS], bool& stop)
+		{
+			ErrorAccumulator e;
+
+			// accumulate squared error
+			for ( double i=-2.5;i<=2.5;i+=0.01 )
+			{
+				// test
+				double actual   = CalculatePiecewise(i,NUM_COEFFICIENTS,testcase);
+				double expected = TargetFunction(i);
+
+				// accumulate error
+				e.AddTestCase(expected,actual);
+			}
+
+			// set stop flag if the error is really low
+			if ( e.GetRMSE() < 0.0000001 )
+			{
+				stop = true;
+			}
+
+			// done
+			return e.GetRMSE();
+		}
+};
+
+
+// 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;
+}
@@ -39,37 +39,29 @@ class PolySearch: public DE::Engine<NUM_COEFFICIENTS>
 {
 	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 TestFitness(const double testcase[NUM_COEFFICIENTS], bool& stop)
 		{
-			double       error      = 0.0;
-			unsigned int test_count = 0;
+			ErrorAccumulator e;

 			// accumulate squared error
 			for ( double i=-1.25;i<=1.25;i+=0.001 )
 			{
 				// test
-				double test   = CalculatePolynomial(i,NUM_COEFFICIENTS,testcase);
-				double actual = TargetFunction(i);
+				double actual   = CalculatePolynomial(i,NUM_COEFFICIENTS,testcase);
+				double expected = TargetFunction(i);

-				// accumulate error^2
-				double e = test - actual;
-				error += e * e;
-
-				// inc test count (divisor to get from SE to MSE)
-				test_count++;
+				// accumulate error
+				e.AddTestCase(expected, actual);
 			}

-			// convert SE to RMSE
-			error = sqrt(error / (double)test_count);
-
 			// set stop flag if the error is really low
-			if ( error < 0.0000001 )
+			if ( e.GetRMSE() < 0.0000001 )
 			{
 				stop = true;
 			}

 			// done
-			return error;
+			return e.GetRMSE();
 		}
 };

@@ -0,0 +1,513 @@
+#include <cstdio>
+#include <cmath>
+#include <cstdlib>
+#include <cassert>
+//#include <vector>
+//#include <algorithm>
+
+
+const double PI = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117;
+
+
+#include <conio.h>
+
+#include "de.h"
+
+
+class PolyDEother: public DE::Engine<5>
+{
+public:
+
+	PolyDEother()
+	{
+		SetRange(-20000000.0,200000000.0);
+	}
+	
+	virtual double CalculateError(const double testcase[5], bool &stop)
+	{
+		double err = 0.0;
+		int    n   = 0;
+
+		for (double x=-10.0;x<=10.0;x+=0.025)
+		{
+			double diff = actual(x) - test(x,testcase);
+			err += (diff * diff);
+			n++;
+		}
+
+		err /= (double)n;
+
+		if ( err < 0.00000000000000000001 )
+		{
+			stop = true;
+		}
+
+		return err;
+	}
+
+	static inline double actual(double x)
+	{
+		const double coeffs[] = { 1.0, 874.0, -451.0, 0.00000000001, -20000000000.0 };
+
+		return test(x,coeffs);
+	}
+
+	static inline double test(double x, const double coeffs[5])
+	{
+		return coeffs[0] + coeffs[1]*x + coeffs[2]*x*x + coeffs[3]*x*x*x + coeffs[4]*x*x*x*x;
+	}
+
+
+
+};
+
+
+#define  DIM  (10)
+
+class PolyDE: public DE::Engine<DIM>
+{
+public:
+
+	double actual_scale0[401*401];
+	double actual_scale1[401*401];
+
+	PolyDE()
+	{
+		SetRange(-0.2,0.2);
+
+		int i=0;
+		for (double t=0.0; t<=1.0; t+=0.0025)
+		{
+			for (double c=0.0; c<=1.0; c+=0.0025)
+			{
+				slerp_actual_t(t,c,actual_scale0[i],actual_scale1[i]);
+				i++;
+				assert(i<(402*402));
+			}
+		}
+
+	}
+
+	virtual double CalculateError(const double testcase[DIM], bool &stop)
+	{
+		double       sumerr = 0.0;
+		unsigned int n      = 0;
+		double err;
+
+		int i = 0;
+		for (double t=0.0; t<=1.0; t+=0.0025)
+		{
+			for (double c=0.0; c<=1.0; c+=0.0025)
+			{
+				double test0,test1;
+				slerp_poly_t(t,c,testcase,test0,test1);
+				
+				err = actual_scale0[i] - test0;
+				sumerr += err * err;
+				err = actual_scale1[i] - test1;
+				sumerr += err * err;
+
+				i++;
+				n++;
+			}
+		}
+
+		sumerr /= (double)n;
+
+		if ( sumerr < 0.0000000000000000000001 || _kbhit() )
+		{
+			stop = true;
+		}
+
+		return sumerr;
+	}
+
+
+
+	static inline void slerp_actual_t(const double& t, const double& costheta, double& scale0, double& scale1)
+	{
+		// calculate coefficients
+		if ( costheta < 1.00000 )
+		{
+			double theta    = acos(costheta);
+			double sintheta = sin(theta);
+			       scale0   = sin((1.0 - t) * theta) / sintheta;
+			       scale1   = sin(t * theta) / sintheta;
+		}
+		else
+		{
+			scale0 = 1.0 - t;
+			scale1 = t;
+		}
+	}
+
+	static inline void slerp_poly_t(const double& t, const double& c, const double coeffs[DIM], double& scale0, double& scale1)
+	{
+		scale0 = sum_poly2(t,c,&coeffs[0]);
+		scale1 = sum_poly2(t,c,&coeffs[DIM>>1]);
+	}
+
+
+
+	static inline double sum_poly2(const double& x, const double& y, const double coeffs[5])
+	{
+		const double x2 = x*x;
+		const double y2 = y*y;
+
+		return coeffs[0] + coeffs[1]*x + coeffs[2]*x2 + coeffs[3]*y + coeffs[4]*y2;
+	}
+	
+	static inline double sum_poly3(const double& x, const double& y, const double coeffs[7])
+	{
+		const double x2 = x*x;
+		const double y2 = y*y;
+		const double x3 = x2*x;
+		const double y3 = y2*y;
+
+		return coeffs[0] + coeffs[1]*x + coeffs[2]*x2 + coeffs[3]*x3 + coeffs[4]*y + coeffs[5]*y2 + coeffs[6]*y3;
+	}
+	
+	static inline double sum_poly4(const double& x, const double& y, const double coeffs[9])
+	{
+		const double x2 = x*x;
+		const double y2 = y*y;
+		const double x3 = x2*x;
+		const double y3 = y2*y;
+		const double x4 = x3*x;
+		const double y4 = y3*y;
+
+		return coeffs[0] + coeffs[1]*x + coeffs[2]*x2 + coeffs[3]*x3 + coeffs[4]*x4 + coeffs[5]*y + coeffs[6]*y2 + coeffs[7]*y3 + coeffs[8]*y4;
+	}
+	
+	static inline double bivariate_poly3(const double& x, const double& y, const double coeffs[10])
+	{
+		const double x2 = x*x;
+		const double y2 = y*y;
+		const double x3 = x2*x;
+		const double y3 = y2*y;
+
+		return x3*(coeffs[9]*y3+coeffs[7]*y2)+x2*(coeffs[8]*y3+coeffs[6]*y2+coeffs[4]*y)+x*(coeffs[5]*y2+coeffs[3]*y+coeffs[1])+coeffs[2]*y+coeffs[0];
+	}
+
+	static inline double bivariate_poly2(const double& x, const double& y, const double coeffs[7])
+	{
+		const double x2 = x*x;
+		const double y2 = y*y;
+		return x2*(coeffs[6]*y2 + coeffs[4]*y) + x*(coeffs[5]*y2 + coeffs[3]*y + coeffs[1]) + coeffs[2]*y + coeffs[0];
+	}
+
+	static inline double bivariate_poly1(const double& x, const double& y, const double coeffs[4])
+	{
+		return x*(coeffs[3]*y + coeffs[1]) + coeffs[2]*y + coeffs[0];
+	}
+
+};
+
+int main()
+{
+	PolyDE* de = new PolyDE;
+
+	de->Solve();
+
+	delete de;
+
+	return 0;
+}
+
+
+
+
+
+
+
+
+
+
+#if 0
+
+class quat
+{
+public:
+
+	double tuple[4];
+
+	static const quat ZERO;
+	static const quat IDENTITY;
+
+	quat()
+	{
+		*this = quat::IDENTITY;
+	}
+
+	~quat()
+	{
+	}
+
+	quat(const quat& other)
+	{
+		*this = other;
+	}
+
+	quat& operator=(const quat& other)
+	{
+		tuple[0] = other.tuple[0];
+		tuple[1] = other.tuple[1];
+		tuple[2] = other.tuple[2];
+		tuple[3] = other.tuple[3];
+		return *this;
+	}
+
+	quat(const double& x, const double& y, const double& z, const double& w)
+	{
+		tuple[0] = x;
+		tuple[1] = y;
+		tuple[2] = z;
+		tuple[3] = w;
+	}
+
+	quat(const double axis[3], const double& angle)
+	{
+		double halfangle = angle / 2.0;
+		double sinangle  = sin(halfangle);
+	
+		tuple[0] = axis[0] * sinangle;
+		tuple[1] = axis[1] * sinangle;
+		tuple[2] = axis[2] * sinangle;
+		tuple[3] = cos(halfangle);
+	}
+	
+	double dot(const quat& other) const
+	{
+		return (tuple[0]*other.tuple[0] + tuple[1]*other.tuple[1] + tuple[2]*other.tuple[2] + tuple[3]*other.tuple[3]);
+	}
+
+	quat& operator*= (double scalar)
+	{
+		tuple[0] *= scalar;
+		tuple[1] *= scalar;
+		tuple[2] *= scalar;
+		tuple[3] *= scalar;
+		return *this;
+	}
+
+	quat operator* (double scalar) const
+	{
+		return quat
+		(
+			tuple[0] * scalar,
+			tuple[1] * scalar,
+			tuple[2] * scalar,
+			tuple[3] * scalar
+		);
+		return *this;
+	}
+
+	quat& operator/= (double scalar)
+	{
+		return (*this *= 1.0/scalar);
+	}
+
+	quat operator/ (double scalar) const
+	{
+		return (*this * 1.0/scalar);
+	}
+
+	double length_squared() const
+	{
+		return dot(*this);
+	}
+
+	double length() const
+	{
+		return sqrt(length_squared());
+	}
+
+	quat& normalize()
+	{
+		return (*this /= length());
+	}
+
+	quat& operator-()
+	{
+		tuple[0] = -tuple[0];
+		tuple[1] = -tuple[1];
+		tuple[2] = -tuple[2];
+		tuple[3] = -tuple[3];
+		return *this;
+	}
+
+	double distance(const quat& other)
+	{
+		quat tmp(tuple[0]-other.tuple[0],tuple[1]-other.tuple[1],tuple[2]-other.tuple[2],tuple[3]-other.tuple[3]);
+		return tmp.length();
+	}
+};
+
+void slerp_coeffs(double t, const quat& q0, const quat& q1, double& scale0, double& scale1)
+{
+	quat from(q0), to(q1);
+	double costheta,sintheta,theta;
+
+	// calculate cosine theta
+	costheta = from.dot(to);
+
+	// adjust signs if necessary
+	if ( costheta < 0.0 )
+	{
+		costheta = -costheta;
+		to       = -to;
+	}
+
+	// calculate coefficients
+	if ( costheta < 1.00000 )
+	{
+		theta    = acos(costheta);
+		sintheta = sin(theta);
+		scale0   = sin((1.0 - t) * theta) / sintheta;
+		scale1   = sin(t * theta) / sintheta;
+	}
+	else
+	{
+		scale0 = 1.0 - t;
+		scale1 = t;
+	}
+
+	//printf("t=%f  s0=%0.10f  s1=%0.10f\n",t,scale0,scale1);
+}
+
+quat slerp(double t, const quat& q0, const quat& q1)
+{
+	double scale0,scale1;
+	slerp_coeffs(t,q0,q1,scale0,scale1);
+
+	// return result
+	quat result
+	(
+		scale0 * q0.tuple[0] + scale1 * q1.tuple[0],
+		scale0 * q0.tuple[1] + scale1 * q1.tuple[1],
+		scale0 * q0.tuple[2] + scale1 * q1.tuple[2],
+		scale0 * q0.tuple[3] + scale1 * q1.tuple[3]
+	);
+
+	result.normalize();
+
+	return result;
+}
+
+quat lerp(double t, const quat& q0, const quat& q1)
+{
+	double scale0 = 1.0 - t;
+	double scale1 = t;
+
+	return quat
+	(
+		scale0 * q0.tuple[0] + scale1 * q1.tuple[0],
+		scale0 * q0.tuple[1] + scale1 * q1.tuple[1],
+		scale0 * q0.tuple[2] + scale1 * q1.tuple[2],
+		scale0 * q0.tuple[3] + scale1 * q1.tuple[3]
+	);
+}
+
+quat slerp_via_lerp_and_normalize(double t, const quat& q0, const quat& q1)
+{
+	double scale0,scale1;
+	slerp_coeffs(t,q0,q1,scale0,scale1);
+
+	double new_t = scale1 / (scale0+scale1);
+
+	quat result = lerp(new_t,q0,q1);
+
+	result.normalize();
+
+	return result;
+}
+
+const quat quat::ZERO     = quat(0.0,0.0,0.0,0.0);
+const quat quat::IDENTITY = quat(0.0,0.0,0.0,1.0);
+
+double vec_length(double vec[3])
+{
+	return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
+}
+
+double random_double(double low, double high)
+{
+	double prn = (double)rand() / (double)RAND_MAX;
+
+	return (prn * (high - low)) + low;
+}
+
+quat random_quat()
+{
+	double axis[3];
+	double angle;
+
+	// random vector elements
+	axis[0] = random_double(-1.0,1.0);
+	axis[1] = random_double(-1.0,1.0);
+	axis[2] = random_double(-1.0,1.0);
+	
+	// normalize vector
+	double veclen = vec_length(axis);
+	axis[0] /= veclen;
+	axis[1] /= veclen;
+	axis[2] /= veclen;
+
+	// random angle
+	angle = random_double(-PI*2.0,PI*2.0);
+
+	// construct a quat
+	return quat(axis,angle);
+}
+
+
+
+bool test()
+{
+	const double MAX_ERROR = 0.0000000000001;
+
+	quat a = random_quat();
+	quat b = random_quat();
+
+	for (double t=0.0; t < 1.01; t+= 0.01)
+	{
+		quat c = slerp(t,a,b);
+		quat d = slerp_via_lerp_and_normalize(t,a,b);
+
+		double error = c.distance(d);
+		if ( error > MAX_ERROR )
+		{
+			printf("actual error = %0.34f, max error = %0.34f\n",error,MAX_ERROR);
+			return false;
+		}
+	}
+
+	return true;
+}
+
+/*
+quat slerp_approx(double t, const quat& q0, const quat& q1)
+{
+	quat from(q0),to(q1);
+
+	// calculate cosine theta
+	double costheta = from.dot(to);
+
+	// adjust signs if necessary
+	if ( costheta < 0.0 )
+	{
+		costheta = -costheta;
+		to       = -to;
+	}
+
+	double new_t = poly[0] + poly[1]*t + poly[2]*t*t + poly[3] + poly[4]*costheta + poly[5]*costheta*costheta;
+
+	quat result = lerp(new_t,q0,q1);
+	result.normalize();
+
+	return result;
+}
+*/
+
+
+#endif
@@ -0,0 +1,158 @@
+#include <cstdio>
+#include <vector>
+#include <algorithm>
+
+#include "de.h"
+
+#define MAX(a,b)          ((a)>(b)?(a):(b))
+#define MIN(a,b)          ((a)<(b)?(a):(b))
+#define CLAMP(min,max,x)  (MAX(min,MIN(max,x)))
+
+class Saboteur: public DE::Engine<16>
+{
+	public:
+		double TestCardMix(int players, int saboteur_cards, int digger_cards)
+		{
+			int total_games = 0;
+			int total_saboteurs = 0;
+			int total_diggers = 0;
+			double saboteur_frequency = 0.0;
+
+			// make deck
+			std::vector<int> deck;
+			for (int i=0;i<saboteur_cards;++i)
+			{
+				deck.push_back(1);
+			}
+			for (int i=0;i<digger_cards;++i)
+			{
+				deck.push_back(0);
+			}
+
+			const int NUM_ITER = 5000;
+			for (int i=0;i<NUM_ITER;++i)
+			{
+				// shuffle deck
+				std::random_shuffle(deck.begin(),deck.end());
+
+				// deal first N cards
+				for (int j=0;j<players;++j)
+				{
+					if ( deck[j] == 1 )
+						++total_saboteurs;
+				}
+			}
+
+			return ((double)total_saboteurs/(double)NUM_ITER) / (double)players;
+		}
+
+		double TestIndividualFitness(int players, int saboteur_cards, int digger_cards, bool verbose=false)
+		{
+			double error_bonus = 0.0;
+			if ( saboteur_cards < 1 || saboteur_cards > 4 )
+				error_bonus += 1000.0;
+			if ( digger_cards < 1 || digger_cards > 7 )
+				error_bonus += 1000.0;
+			if ( saboteur_cards >= players )
+				error_bonus += 500.0;
+			if ( (saboteur_cards+digger_cards) < (players+1) )
+				error_bonus += 500.0;
+			if ( (saboteur_cards+digger_cards) > (players+3) )
+				error_bonus += 500.0;
+
+			if ( verbose )
+			{
+				printf("\n %2d players (%d, %d): ",players,saboteur_cards,digger_cards);
+			}
+
+			if ( error_bonus > 0.0 )
+			{
+				if ( verbose )
+				{
+					printf(" ERROR!  (%f)\n",error_bonus);
+				}
+				return error_bonus;
+			}
+			else
+			{
+				double actual = TestCardMix(players,saboteur_cards,digger_cards);
+				if ( verbose )
+				{
+					printf(" %0.5f\n",actual);
+				}
+				return actual;
+			}
+		}
+
+		double TestFitness(const double testcase[16], bool& stop)
+		{
+			ErrorAccumulator e;
+
+			const int num_players[8] = { 3, 4, 5, 6, 7, 8, 9, 10 };
+			const double expected = 0.363;
+			//const double expected = 1.0/3.0; // saboteurs are 1:3
+
+			for (int i=0;i<8;++i)
+			{
+				int i_saboteur_cards = i * 2;
+				int i_digger_cards   = i_saboteur_cards + 1;
+				int players          = num_players[i];
+
+				int saboteur_cards = (int)testcase[i_saboteur_cards];
+				int digger_cards   = (int)testcase[i_digger_cards];
+
+				double actual = TestIndividualFitness(players,saboteur_cards,digger_cards);
+
+				e.AddTestCase(expected,actual);
+			}
+
+			return e.GetRMSE();
+		}
+};
+
+int main(int argc, char** argv)
+{
+	// init
+	Saboteur s;
+	for (int i=0;i<8;++i)
+	{
+		s.SetRange(  i*2, 1.0, 4.0);
+		s.SetRange(i*2+1, 1.0, 7.0);
+	}
+	s.Reset();
+
+	// print original
+	const double original[16] = { 1,3,1,4,2,4,2,5,3,5,3,6,3,7,4,7 };
+	for (int i=0;i<8;++i)
+	{
+		const int num_players[8] = { 3, 4, 5, 6, 7, 8, 9, 10 };
+		int i_saboteur_cards = i * 2;
+		int i_digger_cards   = i_saboteur_cards + 1;
+		int players          = num_players[i];
+
+		int saboteur_cards = (int)original[i_saboteur_cards];
+		int digger_cards   = (int)original[i_digger_cards];
+
+		s.TestIndividualFitness(players,saboteur_cards,digger_cards,true);
+	}
+
+	// run
+	s.Solve(1000);
+
+	// output
+	const double* best = s.GetBest();
+	for (int i=0;i<8;++i)
+	{
+		const int num_players[8] = { 3, 4, 5, 6, 7, 8, 9, 10 };
+		int i_saboteur_cards = i * 2;
+		int i_digger_cards   = i_saboteur_cards + 1;
+		int players          = num_players[i];
+
+		int saboteur_cards = (int)best[i_saboteur_cards];
+		int digger_cards   = (int)best[i_digger_cards];
+
+		s.TestIndividualFitness(players,saboteur_cards,digger_cards,true);
+	}
+
+	return 0;
+}
@@ -0,0 +1,44 @@
+#include <cstdio>
+
+#include "de.h"
+
+// just a helper function used in calculating the polynomails
+double Polynomial(double x, double a, double b, double c, double d)
+{
+	return a*(x*x*x) + b*(x*x) + c*x + d;
+}
+
+// the function we're searching for
+double TargetFunction(double input)
+{
+	return Polynomial(input,17.0,321.0,-2.0,324.0);
+}
+
+// our DE implementation
+class PolySearch: public DE::Engine<4>
+{
+	public:
+		double TestFitness(const double testcase[4], bool& stop)
+		{
+			double error = 0.0;
+			for ( int i=-100;i<=100;i++ )
+			{
+				double test   = Polynomial((double)i,testcase[0],testcase[1],testcase[2],testcase[3]);
+				double actual = TargetFunction((double)i);
+				error += (test-actual)*(test-actual);
+			}
+			if ( error < 0.00000000000001 )
+			{
+				stop = true;
+			}
+			return error;
+		}
+};
+
+// main program
+int main()
+{
+	PolySearch de;
+	de.Solve();
+	return 0;
+}
Author	SHA1	Message	Date
swaits	69c3fe9d0f	add github mirror	2024-01-29 10:42:37 -07:00
Stephen Waits	a9de6fc165	add README.md with description from previous Bitbucket repo	2020-04-24 08:20:32 -07:00
swaits	3f63ae0e2e	add quaternion search example	2011-01-24 14:18:30 -08:00
swaits	0b158fdb65	add license to individual source files. add intro blurb in de.h.	2010-09-30 16:39:40 -07:00
swaits	9f7cc21ddd	add (MIT) license	2010-09-30 16:19:03 -07:00
swaits	6f81ad879d	add another simple example	2010-09-30 15:53:16 -07:00
swaits	bcbb99be6f	create example dir	2010-09-30 15:39:15 -07:00
swaits	b656b251c4	add saboteur optimization problem. minor compile fixes.	2010-09-30 21:36:37 +00:00
swaits	aac0b98ba9	renamed DE::Engine::CalculateError to DE::Engine::TestFitness added ErrorAccumulator to aid in calculating RMSE converted piecewiseapprox.cpp and polyapprox.cpp to use new pure virtual method name and ErrorAccumulator	2007-08-10 23:43:15 +00:00
swaits	6f03c0b106	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	2007-05-25 01:24:26 +00:00