MaCh3  2.8.0
Reference Guide
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
MCMCBase Class Referenceabstract

Base class for MCMC fitting algorithms. More...

#include <Fitters/Algorithms/MCMCBase.h>

Inheritance diagram for MCMCBase:
[legend]
Collaboration diagram for MCMCBase:
[legend]

Public Member Functions

 MCMCBase (Manager *const fitMan)
 Constructor. More...
 
virtual ~MCMCBase ()=default
 Destructor. More...
 
void RunMCMC () final
 Actual implementation of MCMC fitting algorithm. More...
 
void StartFromPreviousFit (const std::string &FitName) final
 Allow to start from previous fit/chain. More...
 
virtual std::string GetFancyName () const final
 This name is to get include some special algorithm extensions. More...
 
void setChainLength (unsigned int L)
 Set how long chain should be. 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...
 
std::string GetName () const
 Get name of class. More...
 

Protected Member Functions

void DoMCMCStep ()
 The full StartStep->DoStep->EndStep chain. More...
 
virtual void ProposeStep ()=0
 Propose a step. More...
 
void PreStepProcess ()
 Actions before step proposal [start stopwatch]. More...
 
void PostStepProcess ()
 Actions after step proposal [end stopwatch, fill tree]. More...
 
virtual void DoStep ()=0
 The MCMC step proposal and acceptance. More...
 
virtual double AcceptanceProbability ()=0
 Step acceptance probability. More...
 
bool IsStepAccepted (const double acc_prob)
 Is step accepted? More...
 
void AcceptStep ()
 Accept a step. More...
 
void AdaptiveStep ()
 Adaptive MCMC step. More...
 
void PrintProgress (const bool StepsPrint=true)
 Print the progress. More...
 
void CheckAcceptanceRates ()
 Check the acceptance rates. 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...
 
std::vector< double > CalculateBinEdges (double lowerlimit, double upperlimit, int n_points) const
 

Protected Attributes

std::unique_ptr< MulticanonicalMCMCHandler > multicanonicalHandler
 multicanonical handler for umbrella sampling
More...
 
bool out_of_bounds
 Do we reject based on hitting boundaries in systs. More...
 
bool accept
 Accept. More...
 
unsigned int chainLength
 number of steps in chain More...
 
bool anneal
 simulated annealing More...
 
double AnnealTemp
 simulated annealing temperature More...
 
bool multicanonical
 multi-canonical method toggle on/off 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...
 
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

Base class for MCMC fitting algorithms.

Inherits from FitterBase and defines the interface for MCMC-based fitting, including chain management and step handling.

Author
Asher Kaboth

Definition at line 8 of file MCMCBase.h.

Constructor & Destructor Documentation

◆ MCMCBase()

MCMCBase::MCMCBase ( Manager *const  fitMan)

Constructor.

Parameters
fitManA pointer to a manager object, which will handle all settings.

Definition at line 6 of file MCMCBase.cpp.

6  : FitterBase(man) {
7 // *************************
8  // Beginning step number
9  stepStart = 0;
10 
11  // Starting parameters should be thrown
12  out_of_bounds = false;
13  chainLength = Get<unsigned>(fitMan->raw()["General"]["MCMC"]["NSteps"], __FILE__, __LINE__);
14  if (chainLength < 10){
15  MACH3LOG_ERROR("MCMC chain length must be at least 10 steps, otherwise this will result in a floating point exception.");
16  throw MaCh3Exception(__FILE__, __LINE__);
17  }
18 
19  AnnealTemp = GetFromManager<double>(fitMan->raw()["General"]["MCMC"]["AnnealTemp"], -999, __FILE__ , __LINE__);
20  if (AnnealTemp < 0)
21  anneal = false;
22  else
23  {
24  MACH3LOG_INFO("Enabling simulated annealing with T = {}", AnnealTemp);
25  anneal = true;
26  }
27 }
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:37
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:35
FitterBase(Manager *const fitMan)
Constructor.
Definition: FitterBase.cpp:17
Manager * fitMan
The manager for configuration handling.
Definition: FitterBase.h:120
unsigned int stepStart
step start, by default 0 if we start from previous chain then it will be different
Definition: FitterBase.h:133
bool anneal
simulated annealing
Definition: MCMCBase.h:79
unsigned int chainLength
number of steps in chain
Definition: MCMCBase.h:76
double AnnealTemp
simulated annealing temperature
Definition: MCMCBase.h:81
bool out_of_bounds
Do we reject based on hitting boundaries in systs.
Definition: MCMCBase.h:70
Custom exception class used throughout MaCh3.
YAML::Node const & raw() const
Return config.
Definition: Manager.h:48

◆ ~MCMCBase()

virtual MCMCBase::~MCMCBase ( )
virtualdefault

Destructor.

Member Function Documentation

◆ AcceptanceProbability()

virtual double MCMCBase::AcceptanceProbability ( )
protectedpure virtual

Step acceptance probability.

Implemented in MR2T2, and DelayedMR2T2.

◆ AcceptStep()

void MCMCBase::AcceptStep ( )
protected

Accept a step.

Definition at line 250 of file MCMCBase.cpp.

250  {
251 // *************************
252  ++accCount;
253  logLCurr = logLProp;
254 
255  // Loop over systematics and accept
256  for (size_t s = 0; s < systematics.size(); ++s)
257  {
258  systematics[s]->AcceptStep();
259  }
260 }
double logLProp
proposed likelihood
Definition: FitterBase.h:127
int accCount
counts accepted steps
Definition: FitterBase.h:131
double logLCurr
current likelihood
Definition: FitterBase.h:125
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:146

◆ AdaptiveStep()

void MCMCBase::AdaptiveStep ( )
protected

Adaptive MCMC step.

Definition at line 220 of file MCMCBase.cpp.

220  {
221 // *************************
222  // Save the Adaptive output
223  for (const auto &syst : systematics)
224  {
225  if (syst->GetDoAdaption()){
226  syst->UpdateAdaptiveCovariance();
227  }
228  }
229 }

◆ CheckAcceptanceRates()

void MCMCBase::CheckAcceptanceRates ( )
protected

Check the acceptance rates.

Definition at line 172 of file MCMCBase.cpp.

172  {
173 // *******************
174  const auto StepEnd = stepStart + chainLength;
175 
176  // Skip check for reallllly early steps
177  if(step - stepStart < std::min(10.0, static_cast<double>(StepEnd)/10)) {
178  return;
179  }
180  // KS: Do not add "throw", this can crash CI or some short dummy chains for testing
181  if(accCount==0 && step - stepStart > StepEnd/10) {
182  MACH3LOG_CRITICAL("No steps were accepted in the MCMC chain after {} steps. Please check your step sizes.", step - stepStart);
183  }
184 
185  double acc_rate = static_cast<double>(accCount) / static_cast<double>(step - stepStart);
186 
187  if (acc_rate < 0.01) {
188  MACH3LOG_WARN("Acceptance rate is very low: {:.4f}. This may indicate that the proposal distribution is not well-tuned. Consider reducing the step size.", acc_rate);
189  } else if (acc_rate > 0.5) {
190  MACH3LOG_WARN("Acceptance rate is very high: {:.4f}. This may indicate that the proposal distribution is too narrow. Consider increasing the step sizes.", acc_rate);
191  }
192 }
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:38
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:36
unsigned int step
current state
Definition: FitterBase.h:123

◆ DoMCMCStep()

void MCMCBase::DoMCMCStep ( )
protected

The full StartStep->DoStep->EndStep chain.

Starts step timer, prints progress

Step proposal, acceptance etc

Tree filling etc.

Definition at line 101 of file MCMCBase.cpp.

101  {
102 // *******************
104  PreStepProcess();
106  DoStep();
108  PostStepProcess();
109 }
void PreStepProcess()
Actions before step proposal [start stopwatch].
Definition: MCMCBase.cpp:112
virtual void DoStep()=0
The MCMC step proposal and acceptance.
void PostStepProcess()
Actions after step proposal [end stopwatch, fill tree].
Definition: MCMCBase.cpp:126

◆ DoStep()

virtual void MCMCBase::DoStep ( )
protectedpure virtual

The MCMC step proposal and acceptance.

Implemented in MR2T2, and DelayedMR2T2.

◆ GetFancyName()

std::string MCMCBase::GetFancyName ( ) const
finalvirtual

This name is to get include some special algorithm extensions.

Reimplemented from FitterBase.

Definition at line 263 of file MCMCBase.cpp.

263  {
264 // *************************
265  auto originalName = GetName();
266  std::string OutName = "";
267  if(multicanonical){
268  OutName += "MultiCanonical_";
269  }
270 
271  for(unsigned int i = 0; i < systematics.size(); ++i){
272  if(systematics[i]->GetDoAdaption()){
273  OutName += "Adaptive_";
274  break;
275  }
276  }
277 
278  OutName += originalName;
279  return OutName;
280 }
std::string GetName() const
Get name of class.
Definition: FitterBase.h:75
bool multicanonical
multi-canonical method toggle on/off
Definition: MCMCBase.h:84

◆ IsStepAccepted()

bool MCMCBase::IsStepAccepted ( const double  acc_prob)
protected

Is step accepted?

Parameters
acc_probused for telling whether step is accepted or not

Definition at line 232 of file MCMCBase.cpp.

232  {
233 // *************************
234  // Get the random number
235  const double fRandom = M3::rand::Uniform();
236  // Do the accept/reject
237  #ifdef MACH3_DEBUG
238  debugFile << " logLProp: " << logLProp << " logLCurr: " << logLCurr << " acc_prob: " << acc_prob << " fRandom: " << fRandom << std::endl;
239  #endif
240 
241  if (fRandom > acc_prob)
242  {
243  // Reject
244  return false;
245  }
246  return true;
247 }
double Uniform(const double min=0.0, const double max=1.0)
Generate a uniform random number in [min, max).
Definition: Random.h:116

◆ PostStepProcess()

void MCMCBase::PostStepProcess ( )
protected

Actions after step proposal [end stopwatch, fill tree].

Definition at line 126 of file MCMCBase.cpp.

126  {
127 // *************************
128  //KS: Some version of ROOT keep spamming about accessing already deleted object which is wrong and not helpful...
129  int originalErrorLevel = gErrorIgnoreLevel;
130  gErrorIgnoreLevel = kFatal;
131 
132  stepClock->Stop();
133  stepTime = stepClock->RealTime();
134 
135  // Write step to output tree
136  outTree->Fill();
137 
138  // Do Adaptive MCMC
139  AdaptiveStep();
140 
141  if (step % auto_save == 0){
142  outTree->AutoSave();
143  }
144  gErrorIgnoreLevel = originalErrorLevel;
145 }
double stepTime
Time of single step.
Definition: FitterBase.h:153
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
Definition: FitterBase.h:151
int auto_save
auto save every N steps
Definition: FitterBase.h:164
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:162
void AdaptiveStep()
Adaptive MCMC step.
Definition: MCMCBase.cpp:220

◆ PreStepProcess()

void MCMCBase::PreStepProcess ( )
protected

Actions before step proposal [start stopwatch].

Definition at line 112 of file MCMCBase.cpp.

112  {
113 // *******************
114  stepClock->Start();
115  out_of_bounds = false;
116 
117  // Print 10 steps in total
118  if ((step - stepStart) % (chainLength / 10) == 0)
119  {
121  PrintProgress();
122  }
123 }
void CheckAcceptanceRates()
Check the acceptance rates.
Definition: MCMCBase.cpp:172
void PrintProgress(const bool StepsPrint=true)
Print the progress.
Definition: MCMCBase.cpp:149

◆ PrintProgress()

void MCMCBase::PrintProgress ( const bool  StepsPrint = true)
protected

Print the progress.

Parameters
StepsPrintwhether to print info about accepted steps and -LogL

Definition at line 149 of file MCMCBase.cpp.

149  {
150 // *******************
151  if(StepsPrint) MACH3LOG_INFO("Step:\t{}/{}, current: {:.2f}, proposed: {:.2f}", step - stepStart, chainLength, logLCurr, logLProp);
152  if(StepsPrint) MACH3LOG_INFO("Accepted/Total steps: {}/{} = {:.2f}", accCount, step - stepStart, static_cast<double>(accCount) / static_cast<double>(step - stepStart));
153 
154  for (size_t i = 0; i < samples.size(); ++i) {
155  samples[i]->PrintRates();
156  }
157 
158  for (ParameterHandlerBase *cov : systematics) {
159  cov->PrintPreFitCurrPropValues();
160  }
161 #ifdef MACH3_DEBUG
162  if (debug)
163  {
164  debugFile << "\n-------------------------------------------------------" << std::endl;
165  debugFile << "Step:\t" << step + 1 << "/" << chainLength << " | current: " << logLCurr << " proposed: " << logLProp << std::endl;
166  }
167 #endif
168 }
std::vector< SampleHandlerInterface * > samples
Sample holder.
Definition: FitterBase.h:141
Base class for handling systematic uncertainty parameters.

◆ ProposeStep()

virtual void MCMCBase::ProposeStep ( )
protectedpure virtual

Propose a step.

Implemented in MR2T2.

◆ RunMCMC()

void MCMCBase::RunMCMC ( )
finalvirtual

Actual implementation of MCMC fitting algorithm.

Initialise the multicanonical handler

Implements FitterBase.

Definition at line 32 of file MCMCBase.cpp.

32  {
33 // *******************
34  // Multicanonical method toggle from yaml config
35  multicanonical = GetFromManager<bool>(fitMan->raw()["General"]["MCMC"]["Multicanonical"]["Enabled"], false, __FILE__, __LINE__);
36  MACH3LOG_INFO("Multicanonical Method: {}", multicanonical);
37 
38  if (multicanonical) {
40  multicanonicalHandler = std::make_unique<MulticanonicalMCMCHandler>();
41 
42  // Initialize the multicanonical handler with the systematics
43  multicanonicalHandler->InitializeMulticanonicalHandlerConfig(fitMan, systematics);
44 #ifdef MACH3_DEBUG
45  // Enable debug output stream for multicanonical handler if debug is enabled
46  multicanonicalHandler->setDebugStream(&debugFile, debug);
47 #endif
48  }
49 
50  // Save the settings into the output file
51  SaveSettings();
52 
53  // Prepare the output branches
54  PrepareOutput();
55 
56  // Remove obsolete memory and make other checks before fit starts
58 
59  // Print Progress before Propose Step
60  PrintProgress(false);
61 
62  // Only propose step if we are running fresh chain. If we are running from previous chain then it's not needed and we can use usual pipeline
63  if(stepStart == 0) {
64  // Reconfigure the samples, systematics and oscillation for first weight
65  // ProposeStep sets logLProp
66  ProposeStep();
67 
68  // Initialise the value of the multicanonical parameter to the centres of the umbrellas
69  if (multicanonical){
70  multicanonicalHandler->InitializeMulticanonicalParams(systematics);
71  }
72  // Set the current logL to the proposed logL for the 0th step
73  // Accept the first step to set logLCurr: this shouldn't affect the MCMC because we ignore the first N steps in burn-in
75  }
76 
77  // Begin MCMC
78  const auto StepEnd = stepStart + chainLength;
79  for (step = stepStart; step < StepEnd; ++step)
80  {
81  DoMCMCStep();
82  }
83  // Save all the MCMC output
84  SaveOutput();
85 
86  // Process MCMC
87  ProcessMCMC();
88 
89  // Save the adaptive MCMC
90  for (const auto &syst : systematics)
91  {
92  if (syst->GetDoAdaption())
93  {
94  auto adaptive_handler = syst->GetAdaptiveHandler();
95  adaptive_handler->SaveAdaptiveToFile(adaptive_handler->GetOutFileName(), syst->GetName(), true);
96  }
97  }
98 }
void ProcessMCMC()
Process MCMC output.
Definition: FitterBase.cpp:408
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:233
void SaveSettings()
Save the settings that the MCMC was run with.
Definition: FitterBase.cpp:79
void PrepareOutput()
Prepare the output file.
Definition: FitterBase.cpp:155
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:225
std::unique_ptr< MulticanonicalMCMCHandler > multicanonicalHandler
multicanonical handler for umbrella sampling
Definition: MCMCBase.h:67
void DoMCMCStep()
The full StartStep->DoStep->EndStep chain.
Definition: MCMCBase.cpp:101
virtual void ProposeStep()=0
Propose a step.

◆ setChainLength()

void MCMCBase::setChainLength ( unsigned int  L)
inline

Set how long chain should be.

Parameters
Lnew chain length

Definition at line 29 of file MCMCBase.h.

29 { chainLength = L; };

◆ StartFromPreviousFit()

void MCMCBase::StartFromPreviousFit ( const std::string &  FitName)
finalvirtual

Allow to start from previous fit/chain.

Parameters
FitNameName of previous chain

Reimplemented from FitterBase.

Definition at line 195 of file MCMCBase.cpp.

195  {
196 // *******************
197  // Use base class
199 
200  // For MCMC we also need to set stepStart
201  TFile *infile = M3::Open(FitName, "READ", __FILE__, __LINE__);
202  TTree *posts = infile->Get<TTree>("posteriors");
203  unsigned int step_val = 0;
204 
205  posts->SetBranchAddress("step", &step_val);
206  posts->GetEntry(posts->GetEntries() - 1);
207 
208  stepStart = step_val;
209  // KS: Also update number of steps if using adaption
210  for (unsigned int i = 0; i < systematics.size(); ++i) {
211  if (systematics[i]->GetDoAdaption()) {
212  systematics[i]->SetNumberOfSteps(stepStart);
213  }
214  }
215  infile->Close();
216  delete infile;
217 }
virtual void StartFromPreviousFit(const std::string &FitName)
Allow to start from previous fit/chain.
Definition: FitterBase.cpp:350
TFile * Open(const std::string &Name, const std::string &Type, const std::string &File, const int Line)
Opens a ROOT file with the given name and mode.

Member Data Documentation

◆ accept

bool MCMCBase::accept
protected

Accept.

Definition at line 73 of file MCMCBase.h.

◆ anneal

bool MCMCBase::anneal
protected

simulated annealing

Definition at line 79 of file MCMCBase.h.

◆ AnnealTemp

double MCMCBase::AnnealTemp
protected

simulated annealing temperature

Definition at line 81 of file MCMCBase.h.

◆ chainLength

unsigned int MCMCBase::chainLength
protected

number of steps in chain

Definition at line 76 of file MCMCBase.h.

◆ multicanonical

bool MCMCBase::multicanonical
protected

multi-canonical method toggle on/off

Definition at line 84 of file MCMCBase.h.

◆ multicanonicalHandler

std::unique_ptr<MulticanonicalMCMCHandler> MCMCBase::multicanonicalHandler
protected

multicanonical handler for umbrella sampling

Definition at line 67 of file MCMCBase.h.

◆ out_of_bounds

bool MCMCBase::out_of_bounds
protected

Do we reject based on hitting boundaries in systs.

Definition at line 70 of file MCMCBase.h.


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