EvoMinimizer.cc

//************************************************************************//
//                                                                        //
//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)                      //
//                 Julian Pychy (julian@ep1.rub.de)                       //
//                 - Ruhr-Universität Bochum                              //
//                                                                        //
//  This file is part of Pawian.                                          //
//                                                                        //
//  Pawian is free software: you can redistribute it and/or modify        //
//  it under the terms of the GNU General Public License as published by  //
//  the Free Software Foundation, either version 3 of the License, or     //
//  (at your option) any later version.                                   //
//                                                                        //
//  Pawian is distributed in the hope that it will be useful,             //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of        //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         //
//  GNU General Public License for more details.                          //
//                                                                        //
//  You should have received a copy of the GNU General Public License     //
//  along with Pawian.  If not, see <http://www.gnu.org/licenses/>.       //
//                                                                        //
//************************************************************************//

#define mysign(x) (( x > 0 ) - ( x < 0 ))

#include <iostream>
#include <fstream>
#include <boost/random.hpp>
#include "PwaUtils/GlobalEnv.hh"
#include "PwaUtils/EvoMinimizer.hh"
#include "PwaUtils/FitParamsBase.hh"
#include "ErrLogger/ErrLogger.hh"

const double EvoMinimizer::DECREASESIGMAFACTOR = 0.9;
const double EvoMinimizer::INCREASESIGMAFACTOR = 1.1;
const double EvoMinimizer::DECREASELOWTHRESH = 0.2;
const double EvoMinimizer::INCREASEHIGHTHRESH = 0.2;
const double EvoMinimizer::LHSPREADEXIT = 0.01;


// Constructor takes AbsFcn to minimze, start parameters upar, population and iteration
// sizes and the output file name suffix
EvoMinimizer::EvoMinimizer(AbsFcn& theAbsFcn, std::shared_ptr<AbsPawianParameters> upar, int population, int iterations) :
  _population(population)
  , _iterations(iterations)
  , _theAbsFcn(&theAbsFcn)
  , _currentBestParams(theAbsFcn.defaultFitValParms())
  , _defaultFitErrParms(theAbsFcn.defaultFitErrParms())
  , _currentResultFileName("currentEvoResult"+GlobalEnv::instance()->outputFileNameSuffix()+".dat")
{
   // Initialize the best global parameters
   _bestParamsGlobal = upar;
}


// Minimization takes place here
std::vector<double> EvoMinimizer::Minimize(){

   double startlh = (*_theAbsFcn)(_bestParamsGlobal->Params());
   double minlh = startlh;
   int numnoimprovement = 0;

   Info << "Start LH = " << startlh << endmsg;

   for(int i = 0; i<_iterations; i++){

      int numbetterlh = 0;
      double maxitlhspread=0;
      double itlh = minlh;
      _iterationParamBackup = std::shared_ptr<AbsPawianParameters>(_bestParamsGlobal->Clone());
      _bestParamsIteration = std::shared_ptr<AbsPawianParameters>(_bestParamsGlobal->Clone());

      for(int j = 0; j<_population; j++){

         // Get iteration start parameters and shuffle them
	_tmpParams = std::shared_ptr<AbsPawianParameters>(_bestParamsGlobal);
         ShuffleParams();

         // Calc likelihood
         double currentlh = (*_theAbsFcn)(_tmpParams->Params());

         // Get information for break condition
         if(fabs(currentlh - minlh) > maxitlhspread){
               maxitlhspread = fabs(currentlh - minlh);
         }

         // Test for new best lh in iteration;
         if(currentlh < itlh){
            itlh = currentlh;
            _bestParamsIteration = std::shared_ptr<AbsPawianParameters>(_tmpParams->Clone());
         }

         // Count number of lh improvements in the iteration
         if(currentlh < minlh){
            numbetterlh++;
        }
      }

      // Break condition
      if(maxitlhspread < LHSPREADEXIT)
         break;

      // If a new minimum has been found, store information
      // and print parameters
      if(numbetterlh > 0){
	_bestParamsGlobal = std::shared_ptr<AbsPawianParameters>(_bestParamsIteration->Clone());
         minlh = itlh;
         numnoimprovement=0;

        GlobalEnv::instance()->fitParamsBase()->getFitParamVal(_bestParamsGlobal->Params(), _currentBestParams);
        std::ofstream theStream(_currentResultFileName.c_str());
        GlobalEnv::instance()->fitParamsBase()->dumpParams(theStream, _currentBestParams, _defaultFitErrParms);
      }
      else{
         numnoimprovement++;
      }

      // Emergency sigma decrease
      if(numnoimprovement > 5){
         AdjustSigma(DECREASESIGMAFACTOR*DECREASESIGMAFACTOR, 0);
      }

      // Adjust errors depending on ratio of improved likelihoods to population size
      std::string logmessage;
      if(((double)numbetterlh / (double)_population) < DECREASELOWTHRESH){
         AdjustSigma(DECREASESIGMAFACTOR, numbetterlh);
         logmessage = "Decreasing errors.";
      }
      else if(((double)numbetterlh / (double)_population) >= INCREASEHIGHTHRESH){
         AdjustSigma(INCREASESIGMAFACTOR, numbetterlh);
         logmessage = "Increasing errors.";
      }

      // Print iteration summary
      Info << "===============================================" << endmsg;
      Info << "Iteration " << i+1 << " / " << _iterations << " done. Best LH " << minlh << endmsg;
      Info << "Likelihood improvements " << numbetterlh << " / " << _population << endmsg;
      Info << "Likelihood spread " << maxitlhspread << endmsg;
      Info << logmessage << endmsg;
      Info << "===============================================" << endmsg;

   } // Iterations

   return _bestParamsGlobal->Params();
}




// Shuffles parameters using gauss distributions
void EvoMinimizer::ShuffleParams(){

   typedef boost::normal_distribution<double> NormalDistribution;
   typedef boost::mt19937 RandomGenerator;
   typedef boost::variate_generator<RandomGenerator&,	NormalDistribution> GaussianGenerator;
   static RandomGenerator rng(static_cast<unsigned> (time(0)));


   for(unsigned int i=0; i<_tmpParams->Params().size(); i++){

      // Don't touch fixed parameters
      if(_tmpParams->IsFixed(i)){
         continue;
      }

      // Decide whether parameter is increased or decreased
      boost::random::uniform_int_distribution<> coin(0,1);
      bool c = coin(rng);

      // Initialize gaussian width as parameter error
      double sigma = _tmpParams->Error(i);

      // If gaussian collides with parameter limit, reduce width
      if(_tmpParams->HasLimits(i)){
         if(c && (2*sigma > (_tmpParams->UpperLimit(i) - _tmpParams->Value(i)))){
            sigma = (_tmpParams->UpperLimit(i) - _tmpParams->Value(i)) / 2.0;
         }
         else if(!c && (2*sigma > (_tmpParams->Value(i) - _tmpParams->LowerLimit(i)))){
            sigma = (_tmpParams->Value(i) - _tmpParams->LowerLimit(i)) / 2.0;
         }
      }

      // Get random number and set new parameter
      NormalDistribution gaussian_dist(0, sigma);
      GaussianGenerator generator(rng, gaussian_dist);
      double val = fabs(generator());

      if(c) _tmpParams->SetValue(i, _tmpParams->Value(i) + val);
      else  _tmpParams->SetValue(i, _tmpParams->Value(i) - val);

      // Check for limits
      if(_tmpParams->HasLimits(i)){
        if(_tmpParams->Value(i) < _tmpParams->LowerLimit(i))
	    _tmpParams->SetValue(i, _tmpParams->LowerLimit(i));
        if(_tmpParams->Value(i) > _tmpParams->UpperLimit(i))
	    _tmpParams->SetValue(i, _tmpParams->UpperLimit(i));
      }
   }

}



// Increase or decrease parameter errors by a factor
void EvoMinimizer::AdjustSigma(double factor, int numimprovements){

   for(unsigned int i=0; i<_bestParamsGlobal->Params().size(); i++){

      // Don't touch fixed parameters
      if(_bestParamsGlobal->IsFixed(i)){
         continue;
      }

      // When a lh improvement was achieved, don't decrease errors of parameters
      // that changed rapidly, also add a bonus to increasements.
      double afactor=factor;
      double pardiffsigmas = fabs(_bestParamsGlobal->Value(i) - _iterationParamBackup->Value(i)) / _bestParamsGlobal->Error(i);

      if(numimprovements > 0 && pardiffsigmas > 1.5){
         if(factor < 1)        afactor = 1.0;
         else if(factor > 1 )  afactor *= afactor;
      }

      // Set new error
      _bestParamsGlobal->SetError(i, _bestParamsGlobal->Error(i) * afactor);
   }

}