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

#include <Fitters/MinuitFit.h>

Inheritance diagram for MinuitFit:

Collaboration diagram for MinuitFit:

Public Member Functions
	MinuitFit (manager *const fitMan)
	Constructor.

virtual	~MinuitFit ()
	Destructor.

void	RunMCMC () override
	Actual implementation of Minuit Fit algorithm.

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

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

virtual	~LikelihoodFit ()
	Destructor.

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

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

virtual void	RunMCMC ()=0
	Implementation of fitting algorithm.

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

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

virtual	~FitterBase ()
	Destructor for the FitterBase class.

void	AddSampleHandler (SampleHandlerBase *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.

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.

virtual void	RunMCMC ()=0
	The specific fitting algorithm implemented in this function depends on the derived class. It could be Markov Chain Monte Carlo (MCMC), MinuitFit, or another algorithm.

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.

void	RunLLHScan ()
	Perform a 1D likelihood scan.

void	GetStepScaleBasedOnLLHScan ()
	LLH scan is good first estimate of step scale.

void	Run2DLLHScan ()
	Perform a 2D likelihood scan.

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

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

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

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

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

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

void	PrepareOutput ()
	Prepare the output file.

void	SaveOutput ()
	Save output and close files.

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

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

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

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

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

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

bool	fMirroring
	Flag telling if mirroring is used or not.

Protected Attributes inherited from FitterBase
manager *	fitMan
	The manager.

unsigned int	step
	current state

double	logLCurr
	current likelihood

double	logLProp
	proposed likelihood

double	accProb
	current acceptance prob

int	accCount
	counts accepted steps

int	stepStart
	step start if restarting

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

std::vector< double >	syst_llh
	systematic llh breakdowns

std::vector< SampleHandlerBase * >	samples
	Sample holder.

unsigned int	TotalNSamples
	Total number of samples used.

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

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

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

double	stepTime
	Time of single step.

std::unique_ptr< TRandom3 >	random
	Random number.

TFile *	outputFile
	Output.

TDirectory *	CovFolder
	Output cov folder.

TDirectory *	SampleFolder
	Output sample folder.

TTree *	outTree
	Output tree with posteriors.

int	auto_save
	auto save every N steps

bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO.

bool	FileSaved
	Checks if file saved not repeat some operations.

bool	SettingsSaved
	Checks if setting saved not repeat some operations.

bool	OutputPrepared
	Checks if output prepared not repeat some operations.

Detailed Description

Implementation of Minuit fitting algorithm [15].

Author: Kamil Skwarczynski

Definition at line 14 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) {
// *******************
  // 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()));
}

◆ ~MinuitFit()

MinuitFit::~MinuitFit ( )

virtual

Destructor.

Definition at line 33 of file MinuitFit.cpp.

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

Member Function Documentation

◆ GetName()

std::string MinuitFit::GetName ( ) const

inlinevirtual

Get name of class.

Reimplemented from LikelihoodFit.

Definition at line 25 of file MinuitFit.h.

25{return "MinuitFit";};

◆ RunMCMC()

void MinuitFit::RunMCMC ( )

overridevirtual

Actual implementation of Minuit Fit algorithm.

Implements LikelihoodFit.

Definition at line 40 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);
 
  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],
        minuit->SetVariable(ParCounter, ((*it)->GetParName(i)), (*it)->GetParInit(i), (*it)->GetDiagonalError(i)/10);
        minuit->SetVariableValue(ParCounter, (*it)->GetParInit(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)
      {
        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