// @(#)root/tmva $Id: SimulatedAnnealing.cxx 20882 2007-11-19 11:31:26Z rdm $   
// Author: Andreas Hoecker, Joerg Stelzer, Helge Voss

/**********************************************************************************
 * Project: TMVA - a Root-integrated toolkit for multivariate data analysis       *
 * Package: TMVA                                                                  *
 * Class  : SimulatedAnnealing                                                    *
 * Web    : http://tmva.sourceforge.net                                           *
 *                                                                                *
 * Description:                                                                   *
 *      Implementation (see header for description)                               *
 *                                                                                *
 * Authors (alphabetical):                                                        *
 *      Andreas Hoecker <Andreas.Hocker@cern.ch> - CERN, Switzerland              *
 *      Joerg Stelzer   <Joerg.Stelzer@cern.ch>  - CERN, Switzerland              *
 *      Helge Voss      <Helge.Voss@cern.ch>     - MPI-K Heidelberg, Germany      *
 *                                                                                *
 * Copyright (c) 2006:                                                            *
 *      CERN, Switzerland                                                         * 
 *      MPI-K Heidelberg, Germany                                                 * 
 *                                                                                *
 * Redistribution and use in source and binary forms, with or without             *
 * modification, are permitted according to the terms listed in LICENSE           *
 * (http://tmva.sourceforge.net/LICENSE)                                          *
 **********************************************************************************/

//_______________________________________________________________________
//                                                                      
// Implementation of Simulated Annealing fitter  
//_______________________________________________________________________

#include "TRandom.h"
#include "TMath.h"
#include "Riostream.h"
#include "TMVA/SimulatedAnnealing.h"
#include "TMVA/Interval.h"
#include "TMVA/IFitterTarget.h"

ClassImp(TMVA::SimulatedAnnealing)

//_______________________________________________________________________
TMVA::SimulatedAnnealing::SimulatedAnnealing( IFitterTarget& target, const std::vector<Interval*>& ranges )
   : fFitterTarget          ( target ),
     fRandom                ( new TRandom(100) ),
     fRanges                ( ranges ),
     fMaxCalls              ( 500000 ),
     fTemperatureGradient   ( 0.3 ),
     fUseAdaptiveTemperature( kFALSE ),
     fInitialTemperature    ( 1000 ),
     fMinTemperature        ( 0 ),
     fEps                   ( 1e-04 ),
     fNFunLoops             ( 25 ),
     fNEps                  ( 4 ), // needs to be at leas 2 !
     fLogger( "SimulatedAnnealing" )
{   
   // constructor
}

//_______________________________________________________________________
TMVA::SimulatedAnnealing::~SimulatedAnnealing()
{
   // destructor
}

//_______________________________________________________________________
Double_t TMVA::SimulatedAnnealing::Minimize( std::vector<Double_t>& parameters )
{
   // minimisation

   // speed up loops
   UInt_t npar = parameters.size();
      
   // sanity check
   if (npar != fRanges.size()) fLogger << kFATAL << "<Minimize> Mismatch in vector lengths: "
                                       << npar << " != " << fRanges.size() << Endl;

   // set values
   Double_t deltaT = fTemperatureGradient;

   // step width of temperature reduction
   Double_t stepWidth = 3*npar;
   if (npar < 20) stepWidth = 50; // don't go below some minimum 

   UInt_t i; // predefine for compatibility with old platforms (eg, SunOS5)
   std::vector<Double_t> diffFCN; 
   for (i = 0; i < (UInt_t)fNEps; ++i) diffFCN.push_back(1e20);

   // result vector
   std::vector<Double_t> bestParameters( parameters );

   // auxiliary vectors
   std::vector<Double_t> xPars( parameters );
   std::vector<Double_t> yPars( parameters );
   
   // initialize adaptive errors with bold guess
   std::vector<Double_t> adaptiveErrors;
   std::vector<Int_t>    nAccepted;
   for (UInt_t k = 0; k < npar; k++ ) {
      adaptiveErrors.push_back( (fRanges[k]->GetMax() - fRanges[k]->GetMin())/10.0 );
      nAccepted.push_back(0);
   }

   // starting point (note that simulated annealing searches for maximum!)
   Double_t retFCN = - fFitterTarget.EstimatorFunction( xPars );
   Int_t    nCalls = 1;
   Double_t maxFCN = retFCN;
   diffFCN[1] = retFCN;

   // set initial temperature
   Double_t temperature = fInitialTemperature;

   //---------------------------------------------------------------------------------
   // perform the optimization until maximum iteration or required accuracy is reached
   Int_t nIteration = 0;

   // required for adaptive annealing
   std::pair<Int_t,Int_t> optPoint;
   std::pair<Int_t,Int_t> optPoint_previous;

   // begin with new temperature
   Bool_t continueWhile = kTRUE;
   while (continueWhile) {

      // iteration counter
      nIteration++;

      if (fUseAdaptiveTemperature) {
         if (nIteration>0 && ( TMath::Abs( ( (float)nIteration/2. ) - nIteration/2 ) < 0.01 ) ) {
            if (optPoint.first          > 0 && optPoint.second          == 0 && 
                optPoint_previous.first > 0 && optPoint_previous.second == 0) {
               if (stepWidth > 20) {
                  deltaT = TMath::Sqrt( deltaT );
                  stepWidth = 0.5*stepWidth;
               }
            }
         }
      }

      // loop over the iterations before temperature reduction:
      for (Int_t m = 0; m < stepWidth; m++) {

         // adaptive annealing steps
         optPoint_previous.first  = optPoint.first;
         optPoint_previous.second = optPoint.second;
         optPoint.first  = 0;
         optPoint.second = 0;

         // loop over the accepted-function evaluations
         for (Int_t j = 0; j < fNFunLoops; j++) {

            // loop over parameters
            for (UInt_t h = 0; h < npar; h++) {

               // randomize parameter h
               yPars[h] = xPars[h] + gRandom->Uniform(-1.0,1.0)*adaptiveErrors[h];

               // retry if randomising has thrown yPars out of its bounds
               while (yPars[h] < fRanges[h]->GetMin() || yPars[h] > fRanges[h]->GetMax()) {
                  yPars[h] = gRandom->Uniform(-2.0,2.0)*adaptiveErrors[h] + parameters[h];
               }

               // recover previous parameter setting 
               if (h >= 1) yPars[h-1] = xPars[h-1];

               // compute estimator for given variable set (again, searches for maximum)
               Double_t retFCNi = - fFitterTarget.EstimatorFunction( yPars ); 

               // too many function evaluations ? --> stop simulated annealing
               ++nCalls;
               if (nCalls >= fMaxCalls) {
                  for (i = 0; i < npar; i++) parameters[i] = bestParameters[i];
                  return -maxFCN;
               }

               // accept new solution if better FCN value
               if (retFCNi > retFCN) {
                  for (i = 0; i < npar; ++i) xPars[i] = yPars[i];
                  retFCN = retFCNi;
                  ++nAccepted[h];

                  // best FCN value so far, record as new best parameter set
                  if (retFCNi > maxFCN) {
                     for (i = 0; i < npar; ++i) bestParameters[i] = yPars[i];
                     
                     if (m <= stepWidth/2) optPoint.first++;
                     else                  optPoint.second++;
                     
                     maxFCN = retFCNi;
                  }
               } 
               else {

                  // original Metropolis et al. scheme to decide whether a solution is 
                  // accepted if FCN is worse than before;
                  // criterion following:
                  //   N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, 
                  //   "Equation of State Calculations by Fast Computing Machines", 
                  //   J. Chem. Phys., 21, 6, 1087-1092, 1953
                  if (gRandom->Uniform(1.0) < this->GetPerturbationProbability( retFCNi, retFCN, 
                                                                                temperature )) {
                     for (i = 1; i < npar; ++i) xPars[i] = yPars[i];
                     retFCN = retFCNi;
                     ++nAccepted[h];
                  } 
               }
            } 
         } 
         // adjust the adaptiveErrors and the temperature
         // adjust adaptiveErrors to accept between 0.4 and 0.6 of the trial 
         // points at the given temperature
         for (i = 0; i < npar; ++i) {
            Double_t ratio = Double_t(nAccepted[i])/Double_t(fNFunLoops);
            // many FCN improvements, ie, far from minimum --> enhance error
            if      (ratio > 0.6) adaptiveErrors[i] *= (2.0 * (ratio - 0.6) / 0.4 + 1.0);
            // few FCN improvements, ie, closer to minimum --> reduce error
            else if (ratio < 0.4) adaptiveErrors[i] /= (2.0 * (0.4 - ratio) / 0.4 + 1.0);
            // else don't touch the error

            // the error shouldn't be larger than the full variable range
            if (adaptiveErrors[i] > fRanges[i]->GetMax() - fRanges[i]->GetMin()) 
               adaptiveErrors[i] = fRanges[i]->GetMax() - fRanges[i]->GetMin();
         }

         // reset number of accepted functions
         for (i = 0; i < npar; ++i) nAccepted[i] = 0;
      } 

      //  terminate simulated annealing if appropriate 
      diffFCN[1] = retFCN;
      if (TMath::Abs(maxFCN - diffFCN[1]) < fEps) 
         for (i = 0; i < (UInt_t)fNEps; ++i) if (TMath::Abs(retFCN - diffFCN[i]) > fEps) continueWhile = kTRUE;

      // more quite criteria
      if (TMath::Abs(maxFCN) < 1e-06   ) continueWhile = kFALSE;
      if (temperature < fMinTemperature) continueWhile = kFALSE;

      if (continueWhile) {
         // continue annealing
         temperature = deltaT * temperature;
         for (i = fNEps-1; i >= 1; --i) diffFCN[i] = diffFCN[i - 1];
         for (i = 0; i < npar; ++i) xPars[i] = bestParameters[i];
         retFCN = maxFCN;
         printf( "new temp: %g --> maxFCN: %f10.10\n" , temperature, maxFCN );
         fLogger << kVERBOSE << "parameters: ";
         for (UInt_t i=0; i<parameters.size(); i++) 
            cout << parameters[i] << " ";
         fLogger << kVERBOSE << Endl << Endl;
         
      }
   } // end of while loop

   // return best parameter set
   for (i = 0; i < npar; i++ ) parameters[i] = bestParameters[i];

   return - maxFCN;
}

//_______________________________________________________________________
Double_t TMVA::SimulatedAnnealing::GetPerturbationProbability( Double_t E, Double_t Eref, 
                                                               Double_t temperature )
{
   // calculates the probability that a perturbation occured
   return (temperature > 0) ? TMath::Exp( (Eref - E)/temperature ) : 0;
}