//************************************************************************//
// //
// 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 "MinFunctions/EvoMinimizer.hh"
#include "PwaUtils/GlobalEnv.hh"
#include "ErrLogger/ErrLogger.hh"
#include "PwaUtils/GlobalEnv.hh"
#include "ConfigParser/ParserBase.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(std::shared_ptr<AbsFcn> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, int population, int iterations) :
AbsPawianMinimizer(theAbsFcnPtr, upar)
,_population(population)
, _iterations(iterations)
, _evoRatioOfModParams(GlobalEnv::instance()->parser()->evoRatioOfModParams())
, _currentResultFileName("currentEvoResult"+GlobalEnv::instance()->outputFileNameSuffix()+".dat")
{
if (_evoRatioOfModParams <= 0. || _evoRatioOfModParams>1.){
Alert << "_evoRatioOfModParams = " << _evoRatioOfModParams << " not possible\n"
<< "value must be set between 0. and 1. !!!!" << endmsg;
exit(1);
}
Info << "population: " << _population
<<"\niterations: " << _iterations
<< "\nratio of parameters to be modified: " << _evoRatioOfModParams << endmsg;
}
// Minimization takes place here
void EvoMinimizer::minimize(){
double startlh = (*_absFcn)(_bestPawianParams->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>(_bestPawianParams->Clone());
_bestParamsIteration = std::shared_ptr<AbsPawianParameters>(_bestPawianParams->Clone());
for(int j = 0; j<_population; j++){
// Get iteration start parameters and shuffle them
_tmpParams = std::shared_ptr<AbsPawianParameters>(_bestPawianParams->Clone());
ShuffleParams();
// Calc likelihood
double currentlh = (*_absFcn)(_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){
_bestPawianParams = std::shared_ptr<AbsPawianParameters>(_bestParamsIteration->Clone());
minlh = itlh;
numnoimprovement=0;
std::ofstream theStream(_currentResultFileName.c_str());
_bestPawianParams->print(theStream);
}
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();
_minimumReached=true;
}
// 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)));
bool acceptNewParams=false;
while(!acceptNewParams){
for(unsigned int i=0; i<_tmpParams->Params().size(); i++){
// Don't touch fixed parameters
if(_tmpParams->IsFixed(i)){
continue;
}
boost::random::uniform_real_distribution<> coinReal(0.,1.);
double c01 = coinReal(rng);
if (c01 > _evoRatioOfModParams) continue; //only 10% of the parameters changed
acceptNewParams=true;
// 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));
}
}
}
}
void EvoMinimizer::printFitResult(double evtWeightSumData){
if(!_minimumReached){
Alert << "minimum has not been reached!!!" << endmsg;
exit(1);
}
Info <<"Final fit params:" << endmsg;
std::vector<double> pVals=_bestPawianParams->Params();
std::vector<double> pErrs=_bestPawianParams->Errors();
for(unsigned int i=0; i < pVals.size(); ++i){
Info << _bestPawianParams->GetName(i) << " : " << pVals.at(i) << "\t" << pErrs.at(i) << endmsg;
}
double finalLh = (*_absFcn)(_bestPawianParams->Params());
Info <<"final result theLh = "<< finalLh << endmsg;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// calculate AIC, BIC criteria and output selected wave contrib
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
unsigned int noOfFreeFitParams=_bestPawianParams->VariableParameters();
double BICcriterion=2.*finalLh+noOfFreeFitParams*log(evtWeightSumData);
double AICcriterion=2.*finalLh+2.*noOfFreeFitParams;
double AICccriterion=AICcriterion+2.*noOfFreeFitParams*(noOfFreeFitParams+1)/(evtWeightSumData-noOfFreeFitParams-1);
Info << "noOfFreeFitParams:\t" <<noOfFreeFitParams;
Info << "evtWeightSumData:\t" <<evtWeightSumData;
Info << "BIC:\t" << BICcriterion << endmsg;
Info << "AIC:\t" << AICcriterion << endmsg;
Info << "AICc:\t" << AICccriterion << endmsg;
}
// Increase or decrease parameter errors by a factor
void EvoMinimizer::AdjustSigma(double factor, int numimprovements){
for(unsigned int i=0; i<_bestPawianParams->Params().size(); i++){
// Don't touch fixed parameters
if(_bestPawianParams->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(_bestPawianParams->Value(i) - _iterationParamBackup->Value(i)) / _bestPawianParams->Error(i);
if(numimprovements > 0 && pardiffsigmas > 1.5){
if(factor < 1) afactor = 1.0;
else if(factor > 1 ) afactor *= afactor;
}
// Set new error
_bestPawianParams->SetError(i, _bestPawianParams->Error(i) * afactor);
}
}