//************************************************************************//
//                                                                        //
//  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/>.       //
//                                                                        //
//************************************************************************//

#include <iostream>
#include <fstream>
#include <boost/random.hpp>

#include "MinFunctions/MinuitMinimizer.hh"
#include "PwaUtils/GlobalEnv.hh"
#include "FitParams/PwaCovMatrix.hh"
#include "FitParams/MnPawianParameters.hh"
#include "FitParams/ParamDepHandler.hh"

#include "ErrLogger/ErrLogger.hh"
#include "PwaUtils/GlobalEnv.hh"
#include "ConfigParser/ParserBase.hh"

#include "Minuit2/FCNBase.h"
#include "Minuit2/FCNGradientBase.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnUserCovariance.h"

template<typename T>
MinuitMinimizer<T>::MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar) :
  AbsPawianMinimizer<T>(theAbsFcnPtr, upar)
  ,_startMnUserParametersPtr(this->_startPawianParams->mnUserParametersPtr())
  ,_startWithCovMat(false)
{
}

template<typename T>
MinuitMinimizer<T>::MinuitMinimizer(std::shared_ptr<AbsFcn<T>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix) :
  AbsPawianMinimizer<T>(theAbsFcnPtr, upar)
  ,_startMnUserParametersPtr(this->_startPawianParams->mnUserParametersPtr())
  ,_startMnCovMatrix(mnCovMatrix)
  ,_startWithCovMat(true)
{
}


// Minimization takes place here
template<typename T>
void MinuitMinimizer<T>::minimize(){
  //  MnMigrad migrad(*_absFcn, _startPawianParams->mnUserParameters());
  //  MnUserParameters startMnUserP(*_startMnUserParametersPtr);
  unsigned int stratLevel=GlobalEnv::instance()->parser()->minuitStrategyLevel();
  std::shared_ptr<MnMigrad> migrad;
  if (_startWithCovMat) migrad=std::shared_ptr<MnMigrad>(new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr, *_startMnCovMatrix));
  else migrad=std::shared_ptr<MnMigrad>(new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr));
  //  (new MnMigrad(*this->_absFcn, *_startMnUserParametersPtr));
  FunctionMinimum* currentFunctionMinimum=0;

  if(stratLevel==1){
    InfoMsg << "start migrad with strategy level " << 1 << endmsg;
    currentFunctionMinimum= new FunctionMinimum((*migrad)(0, GlobalEnv::instance()->parser()->tolerance()));
  }
  else if(stratLevel==2){
    MnMigrad migrad2a(*this->_absFcn, *_startMnUserParametersPtr, MnStrategy(2));
    InfoMsg << "start migrad with strategy level " << 2 << endmsg;
    currentFunctionMinimum = new FunctionMinimum(migrad2a(0, GlobalEnv::instance()->parser()->tolerance()));
  }
  else{
    Alert << "MIGRAD strategy level " << GlobalEnv::instance()->parser()->minuitStrategyLevel() 
	  << " is not supported!!!" << endmsg;
    exit(1);
  }

 //    MnMigrad migrad(*_absFcn, startMnUserP);
  InfoMsg << "start migrad " << endmsg;
  //  FunctionMinimum currentFunctionMinimum = migrad(0, GlobalEnv::instance()->parser()->tolerance());

  if(currentFunctionMinimum->IsValid()){
    //     return funcMin;
    this->_minimumReached=true;
    _mnFunctionMinimumFinalPtr=std::shared_ptr<FunctionMinimum>(new FunctionMinimum(*currentFunctionMinimum));
    //    _bestPawianParams=std::shared_ptr<AbsPawianParameters>(new MnPawianParameters(_mnFunctionMinimumFinalPtr->UserParameters()));
    this->_bestPawianParams->SetAllValues(_mnFunctionMinimumFinalPtr->UserParameters().Params());
    this->_bestPawianParams->SetAllErrors(_mnFunctionMinimumFinalPtr->UserParameters().Errors());
    ParamDepHandler::instance()->ApplyDependencies(this->_bestPawianParams);
    return;
  }

  // Two more tries to get a valid result unsing strategy 2
  for(int j=0; j<2; j++){
    WarningMsg << "FM is invalid, try with strategy = 2." << endmsg;
    
    // Check minimum covariance matrix
    bool badCovarianceDiagonal=false;
    if(currentFunctionMinimum->HasCovariance()){
      badCovarianceDiagonal = !PwaCovMatrix::DiagonalIsValid(currentFunctionMinimum->UserCovariance());
    }
    
    if(badCovarianceDiagonal){
      InfoMsg << "bad covariance diagonal matrix: Using default errors" << endmsg;
      std::shared_ptr<MnUserParameters> newMnUserParams = this->_startPawianParams->mnUserParametersPtr();
      for(unsigned int i=0; i< currentFunctionMinimum->UserParameters().Params().size();i++){
	newMnUserParams->SetValue(i, currentFunctionMinimum->UserParameters().Params().at(i));
      }
      MnMigrad migrad2(*this->_absFcn, *newMnUserParams, MnStrategy(2));
      if(0!=currentFunctionMinimum) delete currentFunctionMinimum;
      currentFunctionMinimum = new FunctionMinimum(migrad2(0, GlobalEnv::instance()->parser()->tolerance()));
    }
    else{
      std::shared_ptr<MnUserParameters> newMnUserParams = this->_startPawianParams->mnUserParametersPtr();
      for(unsigned int i=0; i< currentFunctionMinimum->UserParameters().Params().size();i++){
	newMnUserParams->SetValue(i, currentFunctionMinimum->UserParameters().Params().at(i));
	newMnUserParams->SetError(i, currentFunctionMinimum->UserParameters().Errors().at(i));
      }
      MnMigrad migrad2(*this->_absFcn, *newMnUserParams, MnStrategy(2));
      if(0!=currentFunctionMinimum) delete currentFunctionMinimum;
      currentFunctionMinimum = new FunctionMinimum(migrad2(0, GlobalEnv::instance()->parser()->tolerance()));
    }
    
    if(currentFunctionMinimum->IsValid()){
      break;
    }
  }
  this->_minimumReached=true;
  _mnFunctionMinimumFinalPtr=std::shared_ptr<FunctionMinimum>(new FunctionMinimum(*currentFunctionMinimum));
  //  _bestPawianParams=std::shared_ptr<AbsPawianParameters>(new MnPawianParameters(_mnFunctionMinimumFinalPtr->UserParameters()));
  this->_bestPawianParams->SetAllValues(_mnFunctionMinimumFinalPtr->UserParameters().Params());
  this->_bestPawianParams->SetAllErrors(_mnFunctionMinimumFinalPtr->UserParameters().Errors());
  ParamDepHandler::instance()->ApplyDependencies(this->_bestPawianParams);
  if(0!=currentFunctionMinimum) delete currentFunctionMinimum;  
}

template<typename T>
void MinuitMinimizer<T>::printFitResult(double evtWeightSumData){
  if(!this->_minimumReached){
    Alert << "minimum has not been reached!!!" << endmsg;
    exit(1);
  }

    double theLh = _mnFunctionMinimumFinalPtr->Fval();

    InfoMsg << "\n\n********************** Final fit parameters *************************" << endmsg;
    this->_bestPawianParams->print(std::cout, true);
    InfoMsg << "\n\n**************** Minuit FunctionMinimum information ******************" << endmsg;
    if(_mnFunctionMinimumFinalPtr->IsValid()) {
      InfoMsg << "\n Function minimum is valid." << endmsg;
    } else {
      InfoMsg << "\n WARNING: Function minimum is invalid!" << endmsg;
    }
    if (_mnFunctionMinimumFinalPtr->HasValidCovariance()) {
      InfoMsg << "\n Covariance matrix is valid." << endmsg;
    } else {
      InfoMsg << "\n WARNING: Covariance matrix is invalid!" << endmsg;
    }
    InfoMsg <<"\n Final LH: "<< std::setprecision(10) << theLh << "\n" << endmsg;
    InfoMsg <<" # of function calls: " << _mnFunctionMinimumFinalPtr->NFcn() << endmsg;
    InfoMsg <<" minimum edm: " << std::setprecision(10) << _mnFunctionMinimumFinalPtr->Edm() << endmsg;
    if (!_mnFunctionMinimumFinalPtr->HasValidParameters()) {
      InfoMsg << " hasValidParameters() returned FALSE" << endmsg;
    }
    if (!_mnFunctionMinimumFinalPtr->HasAccurateCovar()) {
      InfoMsg << " hasAccurateCovar() returned FALSE" << endmsg;
    }
    if (!_mnFunctionMinimumFinalPtr->HasPosDefCovar()) {
      InfoMsg << " hasPosDefCovar() returned FALSE" << endmsg;
      if(_mnFunctionMinimumFinalPtr->HasMadePosDefCovar()) {
	InfoMsg << " hasMadePosDefCovar() returned TRUE" << endmsg;
      }
    }
    if (!_mnFunctionMinimumFinalPtr->HasCovariance()) {
      InfoMsg << " hasCovariance() returned FALSE" << endmsg;
    }
    if (_mnFunctionMinimumFinalPtr->HasReachedCallLimit()) {
      InfoMsg << " hasReachedCallLimit() returned TRUE" << endmsg;
    }
    if (_mnFunctionMinimumFinalPtr->IsAboveMaxEdm()) {
      InfoMsg << " isAboveMaxEdm() returned TRUE" << endmsg;
    }
    if (_mnFunctionMinimumFinalPtr->HesseFailed()) {
      InfoMsg << " hesseFailed() returned TRUE\n" << endmsg;
    }


    ///////////////////////////////////////////////////////////////////////////////////////////////
    // calculate AIC, BIC criteria and output selected wave contrib
    ///////////////////////////////////////////////////////////////////////////////////////////////
    unsigned int noOfFreeFitParams=this->_bestPawianParams->VariableParameters();

    double BICcriterion=2.*theLh+noOfFreeFitParams*log(evtWeightSumData);
    double AICcriterion=2.*theLh+2.*noOfFreeFitParams;
    double AICccriterion=AICcriterion+2.*noOfFreeFitParams*(noOfFreeFitParams+1)
      / (evtWeightSumData-noOfFreeFitParams-1);
    InfoMsg << "noOfFreeFitParams:\t" <<noOfFreeFitParams << endmsg;
    InfoMsg << "evtWeightSumData:\t" <<evtWeightSumData << endmsg;
    InfoMsg << "BIC:\t" << BICcriterion << endmsg;
    InfoMsg << "AIC:\t" << AICcriterion << endmsg;
    InfoMsg << "AICc:\t" << AICccriterion << endmsg;
}

template<typename T>
void MinuitMinimizer<T>::dumpFitResult(){

  std::ostringstream finalResultname;
  std::string outputFileNameSuffix= GlobalEnv::instance()->outputFileNameSuffix();
  finalResultname << "finalResult" << outputFileNameSuffix << ".dat";
  
  std::ofstream theStream ( finalResultname.str().c_str() );
  this->_bestPawianParams->print(theStream);

  //dump covariance matrix
  MnUserCovariance theCovMatrix = _mnFunctionMinimumFinalPtr->UserCovariance();
  std::ostringstream serializationFileName;
  serializationFileName << "serializedOutput" << GlobalEnv::instance()->outputFileNameSuffix() << ".dat";
  std::ofstream serializationStream(serializationFileName.str().c_str());
  boost::archive::text_oarchive boostOutputArchive(serializationStream);
  
  if(_mnFunctionMinimumFinalPtr->HasValidCovariance()){
    const PwaCovMatrix thePwaCovMatrix(theCovMatrix, _mnFunctionMinimumFinalPtr->UserParameters());
    boostOutputArchive << thePwaCovMatrix;
  }
}

template MinuitMinimizer<FCNBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar);
template MinuitMinimizer<FCNGradientBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNGradientBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar);

template MinuitMinimizer<FCNBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix);
template MinuitMinimizer<FCNGradientBase>::MinuitMinimizer(std::shared_ptr<AbsFcn<FCNGradientBase>> theAbsFcnPtr, std::shared_ptr<AbsPawianParameters> upar, std::shared_ptr<MnUserCovariance> mnCovMatrix);
