Implementation of delayed rejection for MR2T2. More...

#include <Fitters/DelayedMR2T2.h>

Inheritance diagram for DelayedMR2T2:

Collaboration diagram for DelayedMR2T2:

Public Member Functions
	DelayedMR2T2 (manager *const fitMan)
	Constructor. More...

virtual	~DelayedMR2T2 ()=default
	Destructor. More...

Public Member Functions inherited from MR2T2
	MR2T2 (manager *const manager)
	Constructor. More...

virtual	~MR2T2 ()=default
	Destructor. More...

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

virtual	~MCMCBase ()=default
	Destructor. More...

void	RunMCMC () override
	Actual implementation of MCMC fitting algorithm. More...

void	StartFromPreviousFit (const std::string &FitName) override
	Allow to start from previous fit/chain. 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 (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. 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	GetStepScaleBasedOnLLHScan ()
	LLH scan is good first estimate of step scale. More...

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

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

void	RunSigmaVarFD ()
	Perform a 1D sigma var for all samples. More...

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

Protected Member Functions
double	AcceptanceProbability () override
	Step acceptance probability. More...

void	DoStep () override
	The MCMC step proposal. More...

void	StoreCurrentStep ()
	Store information about the current proposed step. More...

void	PrepareOutput ()
	override to add in rejection information to output More...

void	ResetSystScale ()
	Reset scale after delay process. More...

void	ScaleSystematics (const double scale)
	Set parameter scale to be original_scale* scale. More...

bool	ProbabilisticDelay () const
	if delay has a chance of occurring More...

Protected Member Functions inherited from MR2T2
void	DoStep () override
	The MCMC step proposal and acceptance. More...

void	ProposeStep () override
	Propose a step. More...

double	AcceptanceProbability () override
	Step acceptance probability. More...

Protected Member Functions inherited from MCMCBase
void	DoMCMCStep ()
	The full StartStep->DoStep->EndStep chain. More...

void	PreStepProcess ()
	Actions before step proposal [start stopwatch]. More...

void	PostStepProcess ()
	Actions after step proposal [end stopwatch, fill tree]. 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 ()
	Print the progress. 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	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. 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)
	For comparison with P-Theta we usually have to apply different parameter values then usual 1, 3 sigma. More...

Protected Attributes
std::vector< std::vector< double > >	current_step_vals
	Store information about the current proposed step. More...

std::vector< double >	start_step_scale
	Stores the initial scale for all parameters. More...

double	initial_scale
	Scale for (non-delayed) step is initial_scale * start_step_scale. More...

double	decay_rate
	How much to decrease the step scale each step. More...

int	max_rejections
	How many rejections allowed? More...

double	MinLogLikelihood
	Minimum (negative) log-likelihood of proposed steps. More...

bool	accepted_delayed
	Was the step we just accepted delayed? More...

bool	delay_on_oob_only
	Delay only if we go out of bounds. More...

double	delay_probability
	Can delay with probability instead. More...

Protected Attributes inherited from MCMCBase
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...

Protected Attributes inherited from FitterBase
manager *	fitMan
	The manager. 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< SampleHandlerBase * >	samples
	Sample holder. More...

unsigned int	TotalNSamples
	Total number of samples used. 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 delayed rejection for MR2T2.

Author: Henry Wallace [13]

Quick Introduction

Delayed rejection is a method to increase the efficiency of MCMC. It is useful if your chain is not mixing well or if adaptive MCMC is struggling. The idea is to propose increasingly small steps centered on the last rejected step.

Configuration

To enable delayed rejection, add the following to your config.yaml file:

Set the algorithm as: General::FittingAlgorithm::DelayedMR2T2

Delayed MCMC will work "as is" but there are additional options which are set in General::MCMC

Scaling Options

Parameters

DecayRate	This should be less than 1 and tells your algorithm how quickly to decrease the step size. By default this is 0.1
MaxRejections	How many steps to delayed-reject before moving on. By default this is 1 and increasing is not recommended
InitialScale	The initial scale of your proposed steps. I.e. propose step with InitialScale*scale in parameter handler. This allows for a really large/small first step to be proposed. Testing finds the default value of 1 to be most effective with adaptive MCMC

Delay Triggers

Parameters

DelayOnlyOutBounds	Here delaying only happens if a step is thrown out of bounds
DelayProbability	In order to mitigate the slowness of delayed MCMC, there is a probability of proposing a delayed step. By default this is 1

Note: Delayed rejection MCMC can significantly improve chain mixing but may increase computational overhead per step.

Definition at line 42 of file DelayedMR2T2.h.

Constructor & Destructor Documentation

◆ DelayedMR2T2()

DelayedMR2T2::DelayedMR2T2 ( manager *const fitMan )

Constructor.

Parameters

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

Definition at line 4 of file DelayedMR2T2.cpp.

                                                   : MR2T2(manager) {
 // *************************
     AlgorithmName = "DelayedMR2T2";
     // Step scale for delayed step = step scale prev delay * decay_rate
     decay_rate = GetFromManager<double>(fitMan->raw()["General"]["MCMC"]["DecayRate"], 0.1);
     // Maximum number of rejections
     max_rejections = GetFromManager<int>(fitMan->raw()["General"]["MCMC"]["MaxRejections"], 1);
     // How big is the first step scale. This scales parameter step by initial_scale * <step scale in ParameterHandler>
     initial_scale = GetFromManager<int>(fitMan->raw()["General"]["MCMC"]["InitialScale"], 1);
  
     // On delayed on out of bounds steps. This saves some compute time and possibly regains efficiency
     delay_on_oob_only = GetFromManager<bool>(fitMan->raw()["General"]["MCMC"]["DelayOnlyOutBounds"], false);
  
     // 
     delay_probability = GetFromManager<double>(fitMan->raw()["General"]["MCMC"]["DelayProbability"], 1.0);
  
  
     if(delay_probability > 1.0){
         MACH3LOG_WARN("Probability of delaying steps > 1, setting to 1");
         delay_probability = 1.0;
     }
     if(delay_probability<=0.0){
         MACH3LOG_CRITICAL("Probability of delaying steps <= 0, setting to 0, delay will no longer have an impact!");
         delay_probability = 0.0;
         // Since we're not delaying at all set max_rejections = 0 to save time
         max_rejections = 0;
     }
  
     MACH3LOG_INFO("Using Delayed MCMC with decay rate: {} and {} allowed rejections", decay_rate, max_rejections);
 }

◆ ~DelayedMR2T2()

virtual DelayedMR2T2::~DelayedMR2T2 ( )

virtualdefault

Destructor.

Member Function Documentation

◆ AcceptanceProbability()

double DelayedMR2T2::AcceptanceProbability ( )

overrideprotectedvirtual

Step acceptance probability.

Implements MCMCBase.

Definition at line 80 of file DelayedMR2T2.cpp.

                                            {
 // *************************
     // From https://arxiv.org/pdf/2010.04190 and DRAM Matlab implementation
     const double numerator = std::max(0.0, std::exp(MinLogLikelihood - logLProp) - 1.0);
     const double denominator = std::exp(MinLogLikelihood - logLCurr) - 1.0;
     if (denominator <= 0.0) return 1.0;
  
     // Large enough that the the minus 1 will make NO difference so the max term cancels
     if(std::isinf(numerator) or std::isinf(denominator) ){
         return MR2T2::AcceptanceProbability();
     }
  
     return std::min(numerator / denominator, 1.0);
 }

◆ DoStep()

void DelayedMR2T2::DoStep ( )

overrideprotectedvirtual

The MCMC step proposal.

OOB Check

Probabalistic condition

Implements MCMCBase.

Definition at line 96 of file DelayedMR2T2.cpp.

                           {
 // *************************
     // Set the initial rejection to false
     StoreCurrentStep();
  
     // Apply initial scale
     ScaleSystematics(initial_scale);
  
     bool accepted = false;
     bool delay = false;
     accepted_delayed = false;
  
     MinLogLikelihood = M3::_LARGE_LOGL_;
  
     for(int i=0; i<=max_rejections; ++i)
     {
         ProposeStep();
  
         // Small hack to ensure we can "leapfrog"
         for (auto &syst : systematics)
         {
             syst->AcceptStep();
         }
  
         if(out_of_bounds or logLProp>MinLogLikelihood){
             continue;
         }
  
         if(i==0){
             accProb = MR2T2::AcceptanceProbability();
         }
         else{
             accProb = AcceptanceProbability();
             delay = true;
         }
         if(IsStepAccepted(accProb)){
             AcceptStep();
             accepted = true;
  
             // Keep track of if step is delayed + accepted for storage
             accepted_delayed = delay;
  
             // Update
             break;
         }
         if(not out_of_bounds and delay_on_oob_only){
             // If we are not out of bounds and only delaying on out of bounds steps
             // We can skip the delay
             break;
         }
         
         if(not ProbabilisticDelay()){
             break;
         }
  
  
         // Temporary, means we can "leapfrog"
         ScaleSystematics(decay_rate);
         MinLogLikelihood = logLProp;
     }
  
     // If we reject we need to reset everything [this is probably a tad too slow...
     // especially in PCA...
     if(!accepted){
         for(int i=0; i<static_cast<int>(systematics.size()); ++i){
             for(int j=0; j<static_cast<int>(systematics[i]->GetParCurrVec().size()); ++j){
                 systematics[i]->SetParCurrProp(j, current_step_vals[i][j]);
             }
         }
     }
     ResetSystScale();
 }

◆ PrepareOutput()

void DelayedMR2T2::PrepareOutput ( )

protected

override to add in rejection information to output

Definition at line 54 of file DelayedMR2T2.cpp.

                                  {
 // *************************
     FitterBase::PrepareOutput();
  
     // Store delayed specific settings
     outTree->Branch("DelayedStep", &accepted_delayed, "DelayedStep/B");
 }

◆ ProbabilisticDelay()

bool DelayedMR2T2::ProbabilisticDelay ( ) const

protected

if delay has a chance of occurring

Definition at line 172 of file DelayedMR2T2.cpp.

                                             {
 // *************************
     // We can delay probabilistically
     return (random->Rndm() > delay_probability);
 }

◆ ResetSystScale()

void DelayedMR2T2::ResetSystScale ( )

protected

Reset scale after delay process.

Definition at line 45 of file DelayedMR2T2.cpp.

                                   {
 // *************************
     for (int i = 0; i < static_cast<int>(systematics.size()); ++i)
     {
         systematics[i]->SetStepScale(start_step_scale[i], false);
     }
 }

◆ ScaleSystematics()

void DelayedMR2T2::ScaleSystematics ( const double scale )

protected

Set parameter scale to be original_scale* scale.

Parameters

scale Multiplicative scale factor

Definition at line 36 of file DelayedMR2T2.cpp.

                                                       {
 // *************************
     for (auto &syst : systematics)
     {
         syst->SetStepScale(scale*syst->GetGlobalStepScale(), false);
     }
 }

◆ StoreCurrentStep()

void DelayedMR2T2::StoreCurrentStep ( )

protected

Store information about the current proposed step.

Definition at line 63 of file DelayedMR2T2.cpp.

                                     {
 // *********************************************
     /*
     Stores values from  step
     */
     // Ensure vector is correctly sized
     current_step_vals.resize(systematics.size());
     start_step_scale.resize(systematics.size());
  
     // Fill with old proposed step, need to copy to avoid using the same step twice
     for(int i=0; i<static_cast<int>(systematics.size()); ++i){
         current_step_vals[i] = std::vector<double>(systematics[i]->GetParCurrVec());
         start_step_scale[i] = systematics[i]->GetGlobalStepScale();
     }
 }

Member Data Documentation

◆ accepted_delayed

bool DelayedMR2T2::accepted_delayed

protected

Was the step we just accepted delayed?

Definition at line 91 of file DelayedMR2T2.h.

◆ current_step_vals

std::vector<std::vector<double> > DelayedMR2T2::current_step_vals

protected

Store information about the current proposed step.

Definition at line 75 of file DelayedMR2T2.h.

◆ decay_rate

double DelayedMR2T2::decay_rate

protected

How much to decrease the step scale each step.

Definition at line 83 of file DelayedMR2T2.h.

◆ delay_on_oob_only

bool DelayedMR2T2::delay_on_oob_only

protected

Delay only if we go out of bounds.

Definition at line 94 of file DelayedMR2T2.h.

◆ delay_probability

double DelayedMR2T2::delay_probability

protected

Can delay with probability instead.

Definition at line 96 of file DelayedMR2T2.h.

◆ initial_scale

double DelayedMR2T2::initial_scale

protected

Scale for (non-delayed) step is initial_scale * start_step_scale.

Definition at line 81 of file DelayedMR2T2.h.

◆ max_rejections

int DelayedMR2T2::max_rejections

protected

How many rejections allowed?

Definition at line 85 of file DelayedMR2T2.h.

◆ MinLogLikelihood

double DelayedMR2T2::MinLogLikelihood

protected

Minimum (negative) log-likelihood of proposed steps.

Definition at line 88 of file DelayedMR2T2.h.

◆ start_step_scale

std::vector<double> DelayedMR2T2::start_step_scale

protected

Stores the initial scale for all parameters.

Definition at line 78 of file DelayedMR2T2.h.

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

Fitters/DelayedMR2T2.h
Fitters/DelayedMR2T2.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Quick Introduction

Configuration

Scaling Options

Delay Triggers

Constructor & Destructor Documentation

◆ DelayedMR2T2()

◆ ~DelayedMR2T2()

Member Function Documentation

◆ AcceptanceProbability()

◆ DoStep()

◆ PrepareOutput()

◆ ProbabilisticDelay()

◆ ResetSystScale()

◆ ScaleSystematics()

◆ StoreCurrentStep()

Member Data Documentation

◆ accepted_delayed

◆ current_step_vals

◆ decay_rate

◆ delay_on_oob_only

◆ delay_probability

◆ initial_scale

◆ max_rejections

◆ MinLogLikelihood

◆ start_step_scale