Changeset View
Changeset View
Standalone View
Standalone View
src/LauMinuit.cc
Show All 20 Lines | |||||
Paul Harrison | Paul Harrison | ||||
Thomas Latham | Thomas Latham | ||||
*/ | */ | ||||
/*! \file LauMinuit.cc | /*! \file LauMinuit.cc | ||||
\brief File containing implementation of LauMinuit methods. | \brief File containing implementation of LauMinuit methods. | ||||
*/ | */ | ||||
#include <iostream> | #include "LauMinuit.hh" | ||||
#include <algorithm> | |||||
#include "TMatrixD.h" | |||||
#include "TVirtualFitter.h" | |||||
#include "LauFitObject.hh" | #include "LauFitObject.hh" | ||||
#include "LauFitter.hh" | #include "LauFitter.hh" | ||||
#include "LauMinuit.hh" | |||||
#include "LauParameter.hh" | #include "LauParameter.hh" | ||||
#include "LauParamFixed.hh" | #include "LauParamFixed.hh" | ||||
#include "TMatrixD.h" | |||||
#include "TVirtualFitter.h" | |||||
#include <algorithm> | |||||
#include <array> | |||||
#include <iostream> | |||||
// It's necessary to define an external function that specifies the address of the function | // It's necessary to define an external function that specifies the address of the function | ||||
// that Minuit needs to minimise. Minuit doesn't know about any classes - therefore | // that Minuit needs to minimise. Minuit doesn't know about any classes - therefore | ||||
// use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this). | // use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this). | ||||
// Here, we use TVirtualFitter* fitter instead of gMinuit, defined below. | // Here, we use TVirtualFitter* fitter instead of gMinuit, defined below. | ||||
// Then, within the external function, invoke an object from this class (LauAllModel), | // Then, within the external function, invoke an object from this class (LauAllModel), | ||||
// and use the member functions to access the parameters/variables. | // and use the member functions to access the parameters/variables. | ||||
extern void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); | extern void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); | ||||
ClassImp(LauMinuit) | ClassImp(LauMinuit) | ||||
LauMinuit::LauMinuit( Int_t maxPar ) : LauAbsFitter(), | LauMinuit::LauMinuit( const UInt_t maxPar, const LauOutputLevel verbosity ) : LauAbsFitter(), | ||||
minuit_(0), | maxPar_{maxPar}, | ||||
maxPar_(maxPar), | outputLevel_{verbosity} | ||||
nParams_(0), | |||||
nFreeParams_(0), | |||||
twoStageFit_(kFALSE), | |||||
useAsymmFitErrors_(kFALSE), | |||||
fitStatus_({-1,0.0,0.0}) | |||||
{ | { | ||||
TVirtualFitter::SetDefaultFitter( "Minuit" ); | TVirtualFitter::SetDefaultFitter( "Minuit" ); | ||||
minuit_ = TVirtualFitter::Fitter( 0, maxPar_ ); | minuit_ = TVirtualFitter::Fitter( nullptr, maxPar_ ); | ||||
} | |||||
LauMinuit::~LauMinuit() | // Set the printout level | ||||
{ | std::array<Double_t,1> argL { static_cast<Double_t>(outputLevel_) }; | ||||
minuit_->ExecuteCommand("SET PRINT", argL.data(), argL.size()); | |||||
if ( outputLevel_ == LauOutputLevel::None ) { | |||||
minuit_->ExecuteCommand("SET NOWARNINGS", argL.data(), 0); | |||||
} | |||||
} | } | ||||
void LauMinuit::initialise( LauFitObject* fitObj, const std::vector<LauParameter*>& parameters ) | void LauMinuit::initialise( LauFitObject* fitObj, const std::vector<LauParameter*>& parameters ) | ||||
{ | { | ||||
// Check whether we're going to use asymmetric errors | // Check whether we're going to use asymmetric errors | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
if (useAsymmFitErrors_ == kTRUE) { | if (useAsymmFitErrors_ == kTRUE) { | ||||
std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl; | std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl; | ||||
std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl; | std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl; | ||||
} else { | } else { | ||||
std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl; | std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl; | ||||
} | } | ||||
} | |||||
// Store the parameters | // Store the parameters | ||||
params_ = parameters; | params_ = parameters; | ||||
// Hook the external likelihood function to this LauFitter::fitter() and this class. | // Hook the external likelihood function to this LauFitter::fitter() and this class. | ||||
minuit_->SetFCN( logLikeFun ); | minuit_->SetFCN( logLikeFun ); | ||||
minuit_->SetObjectFit( fitObj ); | minuit_->SetObjectFit( fitObj ); | ||||
// Clear any stored parameters etc... before using | // Clear any stored parameters etc... before using | ||||
minuit_->Clear(); | minuit_->Clear(); | ||||
nParams_ = params_.size(); | nParams_ = params_.size(); | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl; | std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl; | ||||
std::cout << " : Total number of parameters = " << nParams_ << std::endl; | std::cout << " : Total number of parameters = " << nParams_ << std::endl; | ||||
} | |||||
// Define the default relative error | // Define the default relative error | ||||
const Double_t defaultError(0.01); | const Double_t defaultError(0.01); | ||||
// Set-up the parameters | // Set-up the parameters | ||||
for (UInt_t i = 0; i < nParams_; ++i) { | for (UInt_t i = 0; i < nParams_; ++i) { | ||||
TString name = params_[i]->name(); | TString name = params_[i]->name(); | ||||
Double_t initVal = params_[i]->initValue(); | Double_t initVal = params_[i]->initValue(); | ||||
Show All 9 Lines | for (UInt_t i = 0; i < nParams_; ++i) { | ||||
Double_t minVal = params_[i]->minValue(); | Double_t minVal = params_[i]->minValue(); | ||||
Double_t maxVal = params_[i]->maxValue(); | Double_t maxVal = params_[i]->maxValue(); | ||||
Bool_t secondStage = params_[i]->secondStage(); | Bool_t secondStage = params_[i]->secondStage(); | ||||
if (this->twoStageFit() && secondStage == kTRUE) { | if (this->twoStageFit() && secondStage == kTRUE) { | ||||
params_[i]->fixed(kTRUE); | params_[i]->fixed(kTRUE); | ||||
} | } | ||||
Bool_t fixVar = params_[i]->fixed(); | Bool_t fixVar = params_[i]->fixed(); | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; | std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; | ||||
} | |||||
minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal); | minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal); | ||||
// Fix parameter if required | // Fix parameter if required | ||||
if (fixVar == kTRUE) { | if (fixVar == kTRUE) { | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
std::cout << " : Fixing parameter " << i << std::endl; | std::cout << " : Fixing parameter " << i << std::endl; | ||||
} | |||||
minuit_->FixParameter(i); | minuit_->FixParameter(i); | ||||
} | } | ||||
} | } | ||||
LauParamFixed pred; | LauParamFixed pred; | ||||
nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | ||||
// Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors | // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors | ||||
// for maximum likelihood fit. Very important command, otherwise all | // for maximum likelihood fit. Very important command, otherwise all | ||||
// extracted errors will be too big, and pull distributions will be too narrow! | // extracted errors will be too big, and pull distributions will be too narrow! | ||||
// TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L) | // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L) | ||||
Double_t argL[2]; | |||||
argL[0] = 0.5; | std::array<Double_t,1> argL { 0.5 }; | ||||
fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL, 1); | fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL.data(), argL.size()); | ||||
//argL[0] = 0; | //argL[0] = 0; | ||||
//fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL, 1); | //fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL.data(), argL.size()); | ||||
} | } | ||||
LauFitObject* LauMinuit::getFitObject() | LauFitObject* LauMinuit::getFitObject() | ||||
{ | { | ||||
return (minuit_!=0) ? dynamic_cast<LauFitObject*>( minuit_->GetObjectFit() ) : 0; | return (minuit_!=nullptr) ? dynamic_cast<LauFitObject*>( minuit_->GetObjectFit() ) : nullptr; | ||||
} | } | ||||
const LauAbsFitter::FitStatus& LauMinuit::minimise() | const LauAbsFitter::FitStatus& LauMinuit::minimise() | ||||
{ | { | ||||
Double_t arglist[2]; | std::array<Double_t,2> arglist { | ||||
arglist[0] = 1000*nParams_; // maximum iterations | 1000.0*nParams_, // maximum iterations | ||||
arglist[1] = 0.05; // tolerance -> min EDM = 0.001*tolerance (0.05) | 0.05 // tolerance -> min EDM = 0.001*tolerance (0.05) | ||||
fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist, 2); | }; | ||||
fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist.data(), arglist.size()); | |||||
// Dummy variables - need to feed them to the function | // Dummy variables - need to feed them to the function | ||||
// used for getting NLL, EDM and error matrix status | // used for getting NLL, EDM and error matrix status | ||||
Double_t errdef; | Double_t errdef; | ||||
Int_t nvpar, nparx; | Int_t nvpar, nparx; | ||||
if (fitStatus_.status != 0) { | if (fitStatus_.status != 0) { | ||||
if ( outputLevel_ > LauOutputLevel::None ) { | |||||
std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl; | std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl; | ||||
} | |||||
} else { | } else { | ||||
// Check that the error matrix is ok | // Check that the error matrix is ok | ||||
fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl; | std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl; | ||||
} | |||||
// 0= not calculated at all | // 0= not calculated at all | ||||
// 1= approximation only, not accurate | // 1= approximation only, not accurate | ||||
// 2= full matrix, but forced positive-definite | // 2= full matrix, but forced positive-definite | ||||
// 3= full accurate covariance matrix | // 3= full accurate covariance matrix | ||||
// Fit result was OK. Now get the more precise errors. | // Fit result was OK. Now get the more precise errors. | ||||
fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist, 1); | fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist.data(), 1); | ||||
if (fitStatus_.status != 0) { | if (fitStatus_.status != 0) { | ||||
if ( outputLevel_ > LauOutputLevel::None ) { | |||||
std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl; | std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl; | ||||
} | |||||
} else { | } else { | ||||
// Check that the error matrix is ok | // Check that the error matrix is ok | ||||
fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | ||||
if ( outputLevel_ > LauOutputLevel::Quiet ) { | |||||
std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl; | std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl; | ||||
} | |||||
// 0= not calculated at all | // 0= not calculated at all | ||||
// 1= approximation only, not accurate | // 1= approximation only, not accurate | ||||
// 2= full matrix, but forced positive-definite | // 2= full matrix, but forced positive-definite | ||||
// 3= full accurate covariance matrix | // 3= full accurate covariance matrix | ||||
// Symmetric errors and eror matrix were OK. | // Symmetric errors and eror matrix were OK. | ||||
// Get asymmetric errors if asked for. | // Get asymmetric errors if asked for. | ||||
if (useAsymmFitErrors_ == kTRUE) { | if (useAsymmFitErrors_ == kTRUE) { | ||||
LauFitObject* fitObj = this->getFitObject(); | LauFitObject* fitObj = this->getFitObject(); | ||||
fitObj->withinAsymErrorCalc( kTRUE ); | fitObj->withinAsymErrorCalc( kTRUE ); | ||||
fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist, 1); | fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist.data(), 1); | ||||
fitObj->withinAsymErrorCalc( kFALSE ); | fitObj->withinAsymErrorCalc( kFALSE ); | ||||
if (fitStatus_.status != 0) { | if (fitStatus_.status != 0) { | ||||
std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl; | std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl; | ||||
} | } | ||||
} | } | ||||
} | } | ||||
} | } | ||||
// Print results | // Print results | ||||
fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); | ||||
if ( outputLevel_ > LauOutputLevel::None ) { | |||||
std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl; | std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl; | ||||
} | |||||
// 0= not calculated at all | // 0= not calculated at all | ||||
// 1= approximation only, not accurate | // 1= approximation only, not accurate | ||||
// 2= full matrix, but forced positive-definite | // 2= full matrix, but forced positive-definite | ||||
// 3= full accurate covariance matrix | // 3= full accurate covariance matrix | ||||
if ( outputLevel_ > LauOutputLevel::None ) { | |||||
minuit_->PrintResults(3, fitStatus_.NLL); | minuit_->PrintResults(3, fitStatus_.NLL); | ||||
} | |||||
// Retrieve the covariance matrix from the fitter | // Retrieve the covariance matrix from the fitter | ||||
// For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_ | // For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_ | ||||
// but only the elements within the sub-matrix nFreeParams_ x nFreeParams_ have values, | // but only the elements within the sub-matrix nFreeParams_ x nFreeParams_ have values, | ||||
// the "trailing" elements are zero, so we trim them off. | // the "trailing" elements are zero, so we trim them off. | ||||
Double_t* covMatrix = minuit_->GetCovarianceMatrix(); | Double_t* covMatrix = minuit_->GetCovarianceMatrix(); | ||||
covMatrix_.Clear(); | covMatrix_.Clear(); | ||||
covMatrix_.ResizeTo( nParams_, nParams_ ); | covMatrix_.ResizeTo( nParams_, nParams_ ); | ||||
covMatrix_.SetMatrixArray( covMatrix ); | covMatrix_.SetMatrixArray( covMatrix ); | ||||
covMatrix_.ResizeTo( nFreeParams_, nFreeParams_ ); | covMatrix_.ResizeTo( nFreeParams_, nFreeParams_ ); | ||||
return fitStatus_; | return fitStatus_; | ||||
} | } | ||||
void LauMinuit::fixSecondStageParameters() | void LauMinuit::fixSecondStageParameters() | ||||
{ | { | ||||
for (UInt_t i = 0; i < nParams_; ++i) { | for (UInt_t i{0}; i < nParams_; ++i) { | ||||
Bool_t secondStage = params_[i]->secondStage(); | if ( params_[i]->secondStage() ) { | ||||
if (secondStage == kTRUE) { | |||||
params_[i]->fixed(kTRUE); | params_[i]->fixed(kTRUE); | ||||
minuit_->FixParameter(i); | minuit_->FixParameter(i); | ||||
} | } | ||||
} | } | ||||
LauParamFixed pred; | LauParamFixed pred; | ||||
nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | ||||
} | } | ||||
void LauMinuit::releaseSecondStageParameters() | void LauMinuit::releaseSecondStageParameters() | ||||
{ | { | ||||
for (UInt_t i = 0; i < nParams_; ++i) { | for (UInt_t i{0}; i < nParams_; ++i) { | ||||
Bool_t secondStage = params_[i]->secondStage(); | if ( params_[i]->secondStage() ) { | ||||
if (secondStage == kTRUE) { | |||||
params_[i]->fixed(kFALSE); | params_[i]->fixed(kFALSE); | ||||
minuit_->ReleaseParameter(i); | minuit_->ReleaseParameter(i); | ||||
} | } | ||||
} | } | ||||
LauParamFixed pred; | LauParamFixed pred; | ||||
nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); | ||||
} | } | ||||
void LauMinuit::updateParameters() | void LauMinuit::updateParameters() | ||||
{ | { | ||||
for (UInt_t i = 0; i < nParams_; ++i) { | for (UInt_t i{0}; i < nParams_; ++i) { | ||||
// Get the value and errors from MINUIT | // Get the value and errors from MINUIT | ||||
Double_t value = minuit_->GetParameter(i); | Double_t value { minuit_->GetParameter(i) }; | ||||
Double_t error(0.0); | Double_t error{0.0}; | ||||
Double_t negError(0.0); | Double_t negError{0.0}; | ||||
Double_t posError(0.0); | Double_t posError{0.0}; | ||||
Double_t globalcc(0.0); | Double_t globalcc{0.0}; | ||||
minuit_->GetErrors(i, posError, negError, error, globalcc); | minuit_->GetErrors(i, posError, negError, error, globalcc); | ||||
params_[i]->valueAndErrors(value, error, negError, posError); | params_[i]->valueAndErrors(value, error, negError, posError); | ||||
params_[i]->globalCorrelationCoeff(globalcc); | params_[i]->globalCorrelationCoeff(globalcc); | ||||
} | } | ||||
} | } | ||||
// Definition of the fitting function for Minuit | // Definition of the fitting function for Minuit | ||||
void logLikeFun(Int_t& npar, Double_t* /*first_derivatives*/, Double_t& f, Double_t* par, Int_t /*iflag*/) | void logLikeFun(Int_t& npar, Double_t* /*first_derivatives*/, Double_t& f, Double_t* par, Int_t /*iflag*/) | ||||
{ | { | ||||
// Routine that specifies the negative log-likelihood function for the fit. | // Routine that specifies the negative log-likelihood function for the fit. | ||||
// Used by the MINUIT minimising code. | // Used by the MINUIT minimising code. | ||||
LauFitObject* theModel = LauFitter::fitter()->getFitObject(); | LauFitObject* theModel = LauFitter::fitter().getFitObject(); | ||||
// Set the internal parameters for this model using parameters from Minuit (pars): | // Set the internal parameters for this model using parameters from Minuit (pars): | ||||
theModel->setParsFromMinuit( par, npar ); | theModel->setParsFromMinuit( par, npar ); | ||||
// Set the value of f to be the total negative log-likelihood for the data sample. | // Set the value of f to be the total negative log-likelihood for the data sample. | ||||
f = theModel->getTotNegLogLikelihood(); | f = theModel->getTotNegLogLikelihood(); | ||||
} | } | ||||