mcmc Class Reference

Implementation of MR2T2 algorithm. More...

#include <Fitters/mcmc.h>

Inheritance diagram for mcmc:

Collaboration diagram for mcmc:

Public Member Functions
	mcmc (manager *const fitMan)
	Constructor.
virtual	~mcmc ()
	Destructor.
void	RunMCMC () override
	Actual implementation of MCMC fitting algorithm.
void	setChainLength (unsigned int L)
	Set how long chain should be.
void	StartFromPreviousFit (const std::string &FitName) override
	Allow to start from previous fit/chain.
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 Member Functions
void	ProposeStep ()
	Propose a step.
void	CheckStep ()
	Do we accept the step.
void	PrintProgress ()
	Print the progress.

Private Attributes
bool	reject
	Do we reject based on hitting boundaries in systs.
unsigned int	chainLength
	number of steps in chain
bool	anneal
	simulated annealing
double	AnnealTemp
	simulated annealing temperature

Additional Inherited Members
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 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 MR2T2 algorithm.

Author

Asher Kaboth

Definition at line 7 of file mcmc.h.

Constructor & Destructor Documentation

mcmc()

mcmc::mcmc

(

manager *const

fitMan

)

Constructor.

Parameters

fitMan

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

Definition at line 6 of file mcmc.cpp.

                       : FitterBase(man) {
// *************************
  // Beginning step number
  stepStart = 0;
 
  // Starting parameters should be thrown
  reject = false;
  chainLength = Get<unsigned>(fitMan->raw()["General"]["MCMC"]["NSteps"], __FILE__, __LINE__);
 
  AnnealTemp = GetFromManager<double>(fitMan->raw()["General"]["MCMC"]["AnnealTemp"], -999);
  if(AnnealTemp < 0) anneal = false;
  else
  {
    MACH3LOG_INFO("Enabling simulated annealing with T = {}", AnnealTemp);
    anneal = true;
  }
}

~mcmc()

mcmc::~mcmc ( )

virtual

Destructor.

Definition at line 26 of file mcmc.cpp.

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

Member Function Documentation

CheckStep()

void mcmc::CheckStep ( )

inlineprivate

Do we accept the step.

Definition at line 33 of file mcmc.cpp.

                     {
// **********************
  bool accept = false;
 
  // Set the acceptance probability to zero
  accProb = 0.0;
 
  // Calculate acceptance probability
  if (anneal) accProb = std::min(1., std::exp( -(logLProp-logLCurr) / (std::exp(-step/AnnealTemp))));
  else accProb = std::min(1., std::exp(logLCurr-logLProp));
 
  // Get the random number
  const double fRandom = random->Rndm();
 
  // Do the accept/reject
  if (fRandom <= accProb) {
    accept = true;
    ++accCount;
  } else {
    accept = false;
  }
 
  #ifdef DEBUG
  if (debug) debugFile << " logLProp: " << logLProp << " logLCurr: " << logLCurr << " accProb: " << accProb << " fRandom: " << fRandom << std::endl;
  #endif
 
  // Update all the handlers to accept the step
  if (accept && !reject) {
    logLCurr = logLProp;
 
    // Loop over systematics and accept
    for (size_t s = 0; s < systematics.size(); ++s) {
      systematics[s]->AcceptStep();
    }
  }
 
  stepClock->Stop();
  stepTime = stepClock->RealTime();
 
  // Write step to output tree
  outTree->Fill();
}

GetName()

std::string mcmc::GetName ( ) const

inlinevirtual

Get name of class.

Reimplemented from FitterBase.

Definition at line 27 of file mcmc.h.

{return "MCMC";};

PrintProgress()

void mcmc::PrintProgress ( )

inlineprivate

Print the progress.

Definition at line 206 of file mcmc.cpp.

                         {
// *******************
  MACH3LOG_INFO("Step:\t{}/{}, current: {:.2f}, proposed: {:.2f}", step - stepStart, chainLength, logLCurr, logLProp);
  MACH3LOG_INFO("Accepted/Total steps: {}/{} = {:.2f}", accCount, step - stepStart, static_cast<double>(accCount) / static_cast<double>(step - stepStart));
 
  for (ParameterHandlerBase *cov : systematics) {
    cov->PrintNominalCurrProp();
  }
  #ifdef DEBUG
  if (debug) {
    debugFile << "\n-------------------------------------------------------" << std::endl;
    debugFile << "Step:\t" << step + 1 << "/" << chainLength << "  |  current: " << logLCurr << " proposed: " << logLProp << std::endl;
  }
  #endif
}

ProposeStep()

void mcmc::ProposeStep ( )

inlineprivate

Propose a step.

Definition at line 136 of file mcmc.cpp.

                       {
// *******************
  // Initial likelihood
  double llh = 0.0;
 
  // Initiate to false
  reject = false;
 
  // Loop over the systematics and propose the initial step
  for (size_t s = 0; s < systematics.size(); ++s) {
    // Could throw the initial value here to do MCMC stability studies
    // Propose the steps for the systematics
    systematics[s]->ProposeStep();
 
    // Get the likelihood from the systematics
    syst_llh[s] = systematics[s]->GetLikelihood();
    llh += syst_llh[s];
 
    #ifdef DEBUG
    if (debug) debugFile << "LLH after " << systematics[s]->GetName() << " " << llh << std::endl;
    #endif
  }
 
  // Check if we've hit a boundary in the systematics
  // In this case we can save time by not having to reconfigure the simulation
  if (llh >= M3::_LARGE_LOGL_) {
    reject = true;
    #ifdef DEBUG
    if (debug) debugFile << "Rejecting based on boundary" << std::endl;
    #endif
  }
 
  // Only reweight when we have a good parameter configuration
  // This speeds things up considerably because for every bad parameter configuration we don't have to reweight the MC
  if (!reject)
  {
    // Could multi-thread this
    // But since sample reweight is multi-threaded it's probably better to do that
    for (size_t i = 0; i < samples.size(); ++i)
    {
      samples[i]->Reweight();
    }
 
    //DB for atmospheric event by event sample migration, need to fully reweight all samples to allow event passing prior to likelihood evaluation
    for (size_t i = 0; i < samples.size(); ++i) {
      // Get the sample likelihoods and add them
      sample_llh[i] = samples[i]->GetLikelihood();
      llh += sample_llh[i];
      #ifdef DEBUG
      if (debug) debugFile << "LLH after sample " << i << " " << llh << std::endl;
      #endif
    }
 
  // For when we don't have to reweight, set sample to madness
  } else {
    for (size_t i = 0; i < samples.size(); ++i) {
      // Set the sample_llh[i] to be madly high also to signify a step out of bounds
      sample_llh[i] = M3::_LARGE_LOGL_;
      #ifdef DEBUG
      if (debug) debugFile << "LLH after REJECT sample " << i << " " << llh << std::endl;
      #endif
    }
  }
 
  // Save the proposed likelihood (class member)
  logLProp = llh;
}

RunMCMC()

void mcmc::RunMCMC ( )

overridevirtual

Actual implementation of MCMC fitting algorithm.

Implements FitterBase.

Definition at line 78 of file mcmc.cpp.

                   {
// *******************
  // Save the settings into the output file
  SaveSettings();
 
  // Prepare the output branches
  PrepareOutput();
 
  // Remove obsolete memory and make other checks before fit starts
  SanitiseInputs();
 
  // Reconfigure the samples, systematics and oscillation for first weight
  // ProposeStep sets logLProp
  ProposeStep();
  // Set the current logL to the proposed logL for the 0th step
  // Accept the first step to set logLCurr: this shouldn't affect the MCMC because we ignore the first N steps in burn-in
  logLCurr = logLProp;
 
  // Begin MCMC
  const auto StepEnd = stepStart + chainLength;
  for (step = stepStart; step < StepEnd; ++step)
  {
    stepClock->Start();
    // Set the initial rejection to false
    reject = false;
 
    // Print 10 steps in total
    if ( (step-stepStart) % (chainLength/10) == 0) {
      PrintProgress();
    }
 
    // Propose current step variation and save the systematic likelihood that results in this step being taken
    // Updates logLProp
    ProposeStep();
 
    // Does the MCMC accept this step?
    CheckStep();
 
    // Auto save the output
    if (step % auto_save == 0) outTree->AutoSave();
  }
 
  // Save all the MCMC output
  SaveOutput();
 
  //Save the Adaptive output
  for(const auto& syst : systematics){
    if(syst->GetDoAdaption()){
      syst->GetAdaptiveHandler()->SaveAdaptiveToFile(syst->GetAdaptiveHandler()->GetOutFileName(), syst->GetName(), true);
    }
  }
 
  // Process MCMC
  ProcessMCMC();
}

setChainLength()

void mcmc::setChainLength ( unsigned int  L )

inline

Set how long chain should be.

Definition at line 19 of file mcmc.h.

{ chainLength = L; };

StartFromPreviousFit()

void mcmc::StartFromPreviousFit ( const std::string &  FitName )

overridevirtual

Allow to start from previous fit/chain.

Parameters

FitName

Name of previous chain

Todo:

implement some check that number of params matches etc

Reimplemented from FitterBase.

Definition at line 223 of file mcmc.cpp.

                                                        {
// *******************
  // Use base class
  FitterBase::StartFromPreviousFit(FitName);
 
  // For MCMC we also need to set stepStart
  TFile *infile = new TFile(FitName.c_str(), "READ");
  TTree *posts = infile->Get<TTree>("posteriors");
  int step_val = -1;
 
  posts->SetBranchAddress("step",&step_val);
  posts->GetEntry(posts->GetEntries()-1);
 
  stepStart = step_val + 1;
  // KS: Also update number of steps if using adaption
  for(unsigned int i = 0; i < systematics.size(); ++i){
    if(systematics[i]->GetDoAdaption()){
      systematics[i]->SetNumberOfSteps(step_val);
    }
  }
  infile->Close();
  delete infile;
}

Member Data Documentation

anneal

bool mcmc::anneal

private

simulated annealing

Definition at line 44 of file mcmc.h.

AnnealTemp

double mcmc::AnnealTemp

private

simulated annealing temperature

Definition at line 46 of file mcmc.h.

chainLength

unsigned int mcmc::chainLength

private

number of steps in chain

Definition at line 41 of file mcmc.h.

reject

bool mcmc::reject

private

Do we reject based on hitting boundaries in systs.

Definition at line 39 of file mcmc.h.

---

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

• Fitters/mcmc.h
• Fitters/mcmc.cpp