#include "stdafx.h"
#include "optimizerevolution.h"
#include <boost/random.hpp>

using namespace igi_optimization;

optimizerEvolution::optimizerEvolution(void)
   {
   }

optimizerEvolution::optimizerEvolution(double eps, int maxIter, double startSteps)
   {
   setEvolutionOptimizationParameters(eps, maxIter, startSteps);
   }

optimizerEvolution::~optimizerEvolution(void)
   {
   }

void optimizerEvolution::setEvolutionOptimizationParameters(double eps, int maxIter, double startSteps, double allowWorse)
   {
   if (eps<=0)  throw;  // see what to throw here
   _eps = eps;
   if (maxIter<=0)  throw;  // see what to throw here
   _maxiter = maxIter;
   if (startSteps<=0)  throw;  // see what to throw here
   _startSteps = startSteps;
   _allowWorse = allowWorse;
   }

void optimizerEvolution::setStartPos(valarray<double>& startPos)
   {
   _start.resize(startPos.size());
   _start = startPos;
   }

valarray<double>& optimizerEvolution::startPos()
   {
   return _start;
   }


valarray<double> optimizerEvolution::optimize( objectiveFunction* p_of, double (objectiveFunction::*p_eval)(valarray<double>& pos)) 
   {
   // first some variables that are needed
   int n = p_of->dimOfDefSet();
   // better check if (n > -1);
   // if (n < 1) ;   // trow an exception
   long numfeval = 0;
   valarray<double> distribution(n);
   distribution = _startSteps;

   // use copy of start value
   valarray<double> x = _start;        // was once the parameter at the function call
   double y_best = (p_of->*p_eval)(x); numfeval++;
   // Define a normal random number distribution which produces "double"
   // values between 0 and 1 (0 inclusive, 1 exclusive).
   boost::mt19937 rng;                 // produces randomness out of thin air
   valarray<double> x_new;
   valarray<double> x_best = x;

   boost::normal_distribution<> normal_dist(0,1);  // mean:0 dist:10
   boost::variate_generator<boost::mt19937, boost::normal_distribution<> >
       die(rng, normal_dist);             // glues randomness with mapping

   double allowWorse=_allowWorse;     // this could be a bit bigger like 1.1 .. to walk diagonally....
   // try some points
   while((distribution.max() > _eps) && (numfeval<_maxiter))
      {
      // change values in x randomly
      for (int i=0; i<n; i++)
         {
         // we only change one element at a time...
         // this gives us trouble to walk in diagonal valleys.
         // for this we allow slightly worse positions x_new when allowWorse is >1.
         // all this is only done to addapt the distribution value.
         // The ultimate would be a matrix to define a ellispoid for the variation.
         x_new = x;
         double d = distribution[i];
         double random = die()*d;  
         x_new[i] = x_new[i]+random;
         double y_new = (p_of->*p_eval)(x_new); numfeval++; 
         if (y_new < y_best*allowWorse) // allow slightly bader x values
            {
            x = x_new;
            allowWorse = pow(allowWorse,0.99);
            if (y_new < y_best) //only safe better values
               {
               y_best = y_new; 
               distribution[i] *= 2;   // be optimistic only on better results
               x_best = x_new;
               }
            }
         else
            {
            distribution[i] *= 0.99;
            }

         }
      }

   // returns the optimum :-)
   return x_best;
   }


valarray<double> optimizerEvolution::optimize( objectiveFunction* p_of )
   {
   // call with the typical eval function
   return optimize( p_of, &objectiveFunction::evalAtPos);
   }



QString optimizerEvolution::optimizer_name() const
   {
   return QString("Evolution (simple)");
   }