Implementation of Minuit fitting algorithm [20]. More...

#include <Fitters/MinuitFit.h>

Inheritance diagram for MinuitFit:

Collaboration diagram for MinuitFit:

Public Member Functions
	MinuitFit (Manager *const fitMan)
	Constructor. More...

virtual	~MinuitFit ()
	Destructor. More...

void	RunMCMC () final
	Actual implementation of Minuit Fit algorithm. More...

Public Member Functions inherited from LikelihoodFit
	LikelihoodFit (Manager *const fitMan)
	Constructor. More...

virtual	~LikelihoodFit ()
	Destructor. More...

virtual double	CalcChi2 (const double *x)
	Chi2 calculation over all included samples and syst objects. More...

int	GetNPars () const
	Get total number of params, this sums over all covariance objects. More...

Public Member Functions inherited from FitterBase
	FitterBase (Manager *const fitMan)
	Constructor. More...

virtual	~FitterBase ()
	Destructor for the FitterBase class. More...

void	AddSampleHandler (SampleHandlerInterface *sample)
	This function adds a sample PDF object to the analysis framework. The sample PDF object will be utilized in fitting procedures or likelihood scans. More...

void	AddSystObj (ParameterHandlerBase *cov)
	This function adds a Covariance object to the analysis framework. The Covariance object will be utilized in fitting procedures or likelihood scans. More...

void	DragRace (const int NLaps=100)
	Calculates the required time for each sample or covariance object in a drag race simulation. Inspired by Dan's feature. More...

void	RunLLHScan ()
	Perform a 1D likelihood scan. More...

void	RunLLHMap ()
	Perform a general multi-dimensional likelihood scan. More...

void	GetStepScaleBasedOnLLHScan (const std::string &filename="")
	LLH scan is good first estimate of step scale. More...

void	Run2DLLHScan ()
	Perform a 2D likelihood scan. More...

void	RunSigmaVar ()
	Perform a 1D/2D sigma var for all samples. More...

virtual void	StartFromPreviousFit (const std::string &FitName)
	Allow to start from previous fit/chain. More...

std::string	GetName () const
	Get name of class. More...

Private Attributes
std::unique_ptr< ROOT::Math::Minimizer >	minuit
	Pointer to minimizer, which most often is Minuit. More...

Additional Inherited Members
Protected Member Functions inherited from LikelihoodFit
void	PrepareFit ()
	prepare output and perform sanity checks More...

Protected Member Functions inherited from FitterBase
void	ProcessMCMC ()
	Process MCMC output. More...

void	PrepareOutput ()
	Prepare the output file. More...

void	SaveOutput ()
	Save output and close files. More...

void	SanitiseInputs ()
	Remove obsolete memory and make other checks before fit starts. More...

void	SaveSettings ()
	Save the settings that the MCMC was run with. More...

bool	GetScanRange (std::map< std::string, std::vector< double >> &scanRanges) const
	YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow. More...

void	GetParameterScanRange (const ParameterHandlerBase *cov, const int i, double &CentralValue, double &lower, double &upper, const int n_points, const std::string &suffix="") const
	Helper function to get parameter scan range, central value. More...

bool	CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
	KS: Check whether we want to skip parameter using skip vector. More...

void	CustomRange (const std::string &ParName, const double sigma, double &ParamShiftValue) const
	For comparison with other fitting frameworks (like P-Theta) we usually have to apply different parameter values then usual 1, 3 sigma. More...

Protected Attributes inherited from LikelihoodFit
int	NPars
	Number of all parameters from all covariances. More...

int	NParsPCA
	Number of all parameters from all covariances in PCA base. More...

bool	fMirroring
	Flag telling if mirroring is used or not. More...

Protected Attributes inherited from FitterBase
Manager *	fitMan
	The manager for configuration handling. More...

unsigned int	step
	current state More...

double	logLCurr
	current likelihood More...

double	logLProp
	proposed likelihood More...

double	accProb
	current acceptance prob More...

int	accCount
	counts accepted steps More...

unsigned int	stepStart
	step start, by default 0 if we start from previous chain then it will be different More...

std::vector< double >	sample_llh
	store the llh breakdowns More...

std::vector< double >	syst_llh
	systematic llh breakdowns More...

std::vector< SampleHandlerInterface * >	samples
	Sample holder. More...

unsigned int	TotalNSamples
	Total number of samples used, single SampleHandler can store more than one analysis sample! More...

std::vector< ParameterHandlerBase * >	systematics
	Systematic holder. More...

std::unique_ptr< TStopwatch >	clock
	tells global time how long fit took More...

std::unique_ptr< TStopwatch >	stepClock
	tells how long single step/fit iteration took More...

double	stepTime
	Time of single step. More...

std::unique_ptr< TRandom3 >	random
	Random number. More...

TFile *	outputFile
	Output. More...

TDirectory *	CovFolder
	Output cov folder. More...

TDirectory *	SampleFolder
	Output sample folder. More...

TTree *	outTree
	Output tree with posteriors. More...

int	auto_save
	auto save every N steps More...

bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO. More...

bool	FileSaved
	Checks if file saved not repeat some operations. More...

bool	SettingsSaved
	Checks if setting saved not repeat some operations. More...

bool	OutputPrepared
	Checks if output prepared not repeat some operations. More...

std::string	AlgorithmName
	Name of fitting algorithm that is being used. More...

Detailed Description

Implementation of Minuit fitting algorithm [20].

Author: Kamil Skwarczynski

Definition at line 18 of file MinuitFit.h.

Constructor & Destructor Documentation

◆ MinuitFit()

MinuitFit::MinuitFit ( Manager *const fitMan )

Constructor.

Todo:: KS: Make this in future configurable, for more see: https://root.cern.ch/doc/master/classROOT_1_1Math_1_1Minimizer.html

Definition at line 5 of file MinuitFit.cpp.

                                  : LikelihoodFit(man) {
 // *******************
   AlgorithmName = "MinuitFit";
   // Minimizer type: determines the underlying implementation.
   // Available types include:
   //   - "Minuit2" (recommended modern option)
   //   - "Minuit"  (legacy)
   //   - "Fumili"
   //   - "GSLMultiMin" (for gradient-free minimization)
   //   - "GSLMultiFit"
   //   - "GSLSimAn" (Simulated Annealing)
   const std::string MinimizerType = "Minuit2";
   // Minimizer algorithm (specific to the selected type).
   // For Minuit2, the following algorithms are available:
   //   - "Migrad"   : gradient-based minimization (default)
   //   - "Simplex"  : Nelder-Mead simplex method (derivative-free)
   //   - "Combined" : combination of Simplex and Migrad
   //   - "Scan"     : parameter grid scan
   const std::string MinimizerAlgo = "Migrad";
  
   MACH3LOG_INFO("Creating instance of Minimizer with {} and {}", MinimizerType, MinimizerAlgo);
  
   minuit = std::unique_ptr<ROOT::Math::Minimizer>(
     ROOT::Math::Factory::CreateMinimizer(MinimizerType.c_str(), MinimizerAlgo.c_str()));
  
   std::string Features = gROOT->GetConfigFeatures();
   if (Features.find("minuit2_omp") != std::string::npos) {
     MACH3LOG_ERROR("Minuit in ROOT has been compiled with OMP support");
     MACH3LOG_ERROR("Since MaCh3 uses OMP for reweight this could cause terrible clash");
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
 }

◆ ~MinuitFit()

MinuitFit::~MinuitFit ( )

virtual

Destructor.

Definition at line 41 of file MinuitFit.cpp.

                       {
 // *************************
 }

Member Function Documentation

◆ RunMCMC()

void MinuitFit::RunMCMC ( )

finalvirtual

Actual implementation of Minuit Fit algorithm.

Implements FitterBase.

Definition at line 47 of file MinuitFit.cpp.

                         {
 // *******************
   PrepareFit();
  
   // Remove obsolete memory and make other checks before fit starts
   SanitiseInputs();
  
   //KS: For none PCA this will be equal to normal parameters
   const int NparsMinuitFull = NPars;
   const int NparsMinuit = NParsPCA;
  
   //KS: Set SetFunction we will Minimize
   ROOT::Math::Functor fChi2(this, &MinuitFit::CalcChi2, NparsMinuit);
   minuit->SetFunction(fChi2);
  
   //KS: add config or something
   minuit->SetPrintLevel(2);
   minuit->SetTolerance(0.01);
   minuit->SetMaxFunctionCalls(fitMan->raw()["General"]["Minuit2"]["NSteps"].as<unsigned>());
   minuit->SetMaxIterations(10000);
  
  
   // Save the adaptive MCMC
   for (const auto &syst : systematics)
   {
     if (syst->GetDoAdaption())
     {
       MACH3LOG_ERROR("Param Handler {} has enabled Adaption, this is not needed for Minuit so please turn it off", syst->GetDoAdaption());
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
   }
  
   MACH3LOG_INFO("Preparing Minuit");
   int ParCounter = 0;
   for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
   {
     if(!(*it)->IsPCA())
     {
       for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
       {
         // KS: Index, name, prior, step scale [different to MCMC],
         // divide by 10 because this is what BANFF used to do...
         minuit->SetVariable(ParCounter, ((*it)->GetParName(i)), (*it)->GetParPreFit(i), (*it)->GetDiagonalError(i)/10);
         minuit->SetVariableValue(ParCounter, (*it)->GetParPreFit(i));
         //KS: lower bound, upper bound, if Mirroring enabled then ignore
         if(!fMirroring) minuit->SetVariableLimits(ParCounter, (*it)->GetLowerBound(i), (*it)->GetUpperBound(i));
         if((*it)->IsParameterFixed(i))
         {
           minuit->FixVariable(ParCounter);
         }
       }
     }
     else
     {
       for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
       {
         // KS: Index, name, prior, step scale [different to MCMC],
         // divide by 10 because this is what BANFF used to do...
         minuit->SetVariable(ParCounter, Form("%i_PCA", i), (*it)->GetPCAHandler()->GetParPropPCA(i), (*it)->GetPCAHandler()->GetEigenValuesMaster()[i]/10);
         if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
         {
           minuit->FixVariable(ParCounter);
         }
       }
     }
   }
  
   minuit->SetPrintLevel(2);
  
   MACH3LOG_INFO("Starting MIGRAD");
   minuit->Minimize();
  
   MACH3LOG_INFO("Starting HESSE");
   minuit->Hesse();
   outputFile->cd();
   
   TVectorD* MinuitParValue = new TVectorD(NparsMinuitFull);
   TVectorD* MinuitParError = new TVectorD(NparsMinuitFull);
   TMatrixDSym* Postmatrix = new TMatrixDSym(NparsMinuitFull);
  
   for(int i = 0; i < NparsMinuitFull; ++i)
   {
     (*MinuitParValue)(i) = 0;
     (*MinuitParError)(i) = 0;
     for(int j = 0; j < NparsMinuitFull; ++j)
     {
       (*Postmatrix)(i,j) = 0;
       (*Postmatrix)(i,j) = minuit->CovMatrix(i,j);
     }
   }
  
   ParCounter = 0;
   const double *X = minuit->X();
   const double *err = minuit->Errors();
   for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
   {
     if(!(*it)->IsPCA())
     {
       for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
       {
         double ParVal = X[ParCounter];
         //KS: Basically apply mirroring for parameters out of bounds
         if(fMirroring)
         {
           if(ParVal < (*it)->GetLowerBound(i))
           {
             ParVal = (*it)->GetLowerBound(i) + ((*it)->GetLowerBound(i) - ParVal);
           }
           else if (ParVal > (*it)->GetUpperBound(i))
           {
             ParVal = (*it)->GetUpperBound(i) - ( ParVal - (*it)->GetUpperBound(i));
           }
         }
         (*MinuitParValue)(ParCounter) = ParVal;
         (*MinuitParError)(ParCounter) = err[ParCounter];
         //KS: For fixed params HESS will not calculate error so we need to pass prior error
         if((*it)->IsParameterFixed(i))
         {
           (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
           (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
         }
       }
     }
     else
     {
       //KS: We need to convert parameters from PCA to normal base
       TVectorD ParVals((*it)->GetNumParams());
       TVectorD ParVals_PCA((*it)->GetNParameters());
  
       TVectorD ErrorVals((*it)->GetNumParams());
       TVectorD ErrorVals_PCA((*it)->GetNParameters());
  
       TMatrixD MatrixVals((*it)->GetNumParams(), (*it)->GetNumParams());
       TMatrixD MatrixVals_PCA((*it)->GetNParameters(), (*it)->GetNParameters());
  
       //First save them
       //KS: This code is super convoluted as MaCh3 can store separate matrices while Minuit has one matrix. In future this will be simplified, keep it like this for now.
       const int StartVal = ParCounter;
       for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
       {
         ParVals_PCA(i) = X[ParCounter];
         ErrorVals_PCA(i) = err[ParCounter];
         int ParCounterMatrix = StartVal;
         for(int j = 0; j < (*it)->GetNParameters(); ++j, ++ParCounterMatrix)
         {
           MatrixVals_PCA(i,j) = minuit->CovMatrix(ParCounter,ParCounterMatrix);
         }
       }
       ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
       ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
       MatrixVals.Mult(((*it)->GetPCAHandler()->GetTransferMatrix()),MatrixVals_PCA);
  
       ParCounter = StartVal;
       //KS: Now after going from PCA to normal let';s save it
       for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
       {
         (*MinuitParValue)(ParCounter) = ParVals(i);
         (*MinuitParError)(ParCounter) = std::fabs(ErrorVals(i));
         int ParCounterMatrix = StartVal;
         for(int j = 0; j < (*it)->GetNumParams(); ++j, ++ParCounterMatrix)
         {
           (*Postmatrix)(ParCounter,ParCounterMatrix)  = MatrixVals(i,j);
         }
         //If fixed take prior
         if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
         {
           (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
           (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
         }
       }
     }
   }
  
   MinuitParValue->Write("MinuitParValue");
   MinuitParError->Write("MinuitParError");
   Postmatrix->Write("Postmatrix");
   delete MinuitParValue;
   delete MinuitParError;
   delete Postmatrix;
   // Save all the output
   SaveOutput();
 }

Member Data Documentation

◆ minuit

std::unique_ptr<ROOT::Math::Minimizer> MinuitFit::minuit

private

Pointer to minimizer, which most often is Minuit.

Definition at line 29 of file MinuitFit.h.

The documentation for this class was generated from the following files:

Fitters/MinuitFit.h
Fitters/MinuitFit.cpp

Public Member Functions

Private Attributes

Additional Inherited Members