Contains information about adaptive covariance matrix [15] [10]. More...

#include <Parameters/AdaptiveMCMCHandler.h>

Collaboration diagram for AdaptiveMCMCHandler:

Public Member Functions
	AdaptiveMCMCHandler ()
	Constructor. More...

virtual	~AdaptiveMCMCHandler ()
	Destructor. More...

void	Print () const
	Print all class members. More...

bool	InitFromConfig (const YAML::Node &adapt_manager, const std::string &matrix_name_str, const std::vector< std::string > parameter_names, const std::vector< double > parameters, const std::vector< double > fixed, const std::vector< bool > param_skip_adapt_flags, const TMatrixDSym *throwMatrix, const double initial_scale_)
	Read initial values from config file. More...

void	CreateNewAdaptiveCovariance ()
	If we don't have a covariance matrix to start from for adaptive tune we need to make one! More...

void	SetAdaptiveBlocks (const std::vector< std::vector< int >> &block_indices)
	HW: sets adaptive block matrix. More...

void	SaveAdaptiveToFile (const std::string &outFileName, const std::string &systematicName, const bool is_final=false)
	HW: Save adaptive throw matrix to file. More...

void	SetThrowMatrixFromFile (const std::string &matrix_file_name, const std::string &matrix_name, const std::string &means_name, bool &use_adaptive)
	sets throw matrix from a file More...

void	CheckMatrixValidityForAdaption (const TMatrixDSym *covMatrix) const
	Check if there are structures in matrix that could result in failure of adaption fits. this mostly cover cases of highly correlated block which seem to results in adaption failure. More...

void	UpdateAdaptiveCovariance ()
	Method to update adaptive MCMC [15]. More...

bool	IndivStepScaleAdapt () const
	Tell whether we want reset step scale or not. More...

bool	UpdateMatrixAdapt ()
	Tell whether matrix should be updated. More...

bool	AdaptionUpdate () const
	To be fair not a clue... More...

bool	SkipAdaption () const
	Tell if we are Skipping Adaption. More...

int	GetParIndex (const std::string &name) const
	Get the index of the parameter given its name. More...

std::vector< int >	GetParIndicesFromNames (const std::vector< std::string > &param_names) const
	Convert vector of parameter names to indices. More...

void	SetParams (const std::vector< double > *params)
	Set the current values of the parameters. More...

void	SetFixed (const std::vector< double > *fix)
	Set the fixed parameters. More...

int	GetNFixed () const
	Get the number of fixed parameters. More...

int	GetNumParams () const
	Get the current values of the parameters. More...

bool	IsFixed (const int ipar) const
	Check if a parameter is fixed. More...

double	CurrVal (const int par_index) const
	Get Current value of parameter. More...

int	GetTotalSteps () const
	Get Total Number of Steps. More...

void	SetTotalSteps (const int nsteps)
	Change Total Number of Steps to new value. More...

void	IncrementNSteps ()
	Increase by one number of total steps. More...

void	IncrementAcceptedSteps ()

const TMatrixDSym *	GetAdaptiveCovariance () const
	Return matrix used for adaption. More...

std::vector< double >	GetParameterMeans () const
	Get the parameter means used in the adaptive handler. More...

std::string	GetOutFileName () const
	Get Name of Output File. More...

double	GetAdaptionScale () const
	Get the current adaption scale. More...

bool	GetUseRobbinsMonro () const
	Use Robbins-Monro approach? More...

void	UpdateRobbinsMonroScale ()
	Update the scale factor for Robbins-Monro adaption. More...

Private Member Functions
void	CalculateRobbinsMonroStepLength ()
	Calculate the constant step length for Robbins-Monro adaption. More...

Private Attributes
int	start_adaptive_throw

int	start_adaptive_update
	When do we stop update the adaptive matrix. More...

int	end_adaptive_update
	Steps between changing throw matrix. More...

int	adaptive_update_step
	Steps between changing throw matrix. More...

int	adaptive_save_n_iterations

std::string	output_file_name
	Name of the file to save the adaptive matrices into. More...

std::vector< int >	adapt_block_matrix_indices
	Indices for block-matrix adaption. More...

std::vector< int >	adapt_block_sizes
	Size of blocks for adaption. More...

std::vector< double >	par_means
	Mean values for all parameters. More...

TMatrixDSym *	adaptive_covariance
	Full adaptive covariance matrix. More...

int	total_steps
	Total number of MCMC steps. More...

double	adaption_scale
	Scaling factor. More...

double	initial_scale
	Initial scaling factor. More...

const std::vector< double > *	_fFixedPars
	Vector of fixed parameters. More...

const std::vector< double > *	_fCurrVal
	Current values of parameters. More...

const std::vector< std::string > *	_ParamNames
	Parameter names. More...

int	acceptance_rate_batch_size
	Acceptance rate in the current batch. More...

bool	use_robbins_monro
	Use Robbins Monro https://arxiv.org/pdf/1006.3690. More...

double	target_acceptance
	Target acceptance rate for Robbins Monro. More...

double	c_robbins_monro
	Constant "step scaling" factor for Robbins-Monro. More...

int	n_rm_restarts
	Number of restarts for Robbins Monro (so far) More...

int	total_rm_restarts
	Total number of restarts ALLOWED for Robbins Monro. More...

bool	prev_step_accepted
	Need to keep track of whether previous step was accepted for RM. More...

const std::vector< bool > *	_param_skip_adapt_flags
	Parameters to skip during adaption. More...

const TMatrixDSym *	initial_throw_matrix
	Initial throw matrix. More...

bool	initial_throw_matrix_saved = false
	Flag to check if initial throw matrix has been saved. More...

Detailed Description

Contains information about adaptive covariance matrix [15] [10].

Author: Henry Wallace; Hank Hua

Adaptation can be performed in two ways:

Haario/Roberts-Rosenthal style: Uses a constant scale factor of
\[ \frac{2.38}{\sqrt{\text{num\_params}}} \]
Robbins-Monro style: Dynamically updates the scale factor to target a specific acceptance rate, potentially improving mixing.

Definition at line 25 of file AdaptiveMCMCHandler.h.

Constructor & Destructor Documentation

◆ AdaptiveMCMCHandler()

AdaptiveMCMCHandler::AdaptiveMCMCHandler ( )

Constructor.

Robbins-Monro adaption

Definition at line 4 of file AdaptiveMCMCHandler.cpp.

                                          {
 // ********************************************
   start_adaptive_throw  = 0;
   start_adaptive_update = 0;
   end_adaptive_update   = 1;
   adaptive_update_step  = 1000;
   total_steps = 0;
   par_means = {};
   adaptive_covariance = nullptr;
  
   use_robbins_monro = false;
   total_rm_restarts = 0;
   
   target_acceptance = 0.234;
   acceptance_rate_batch_size = 10000;
 }

◆ ~AdaptiveMCMCHandler()

AdaptiveMCMCHandler::~AdaptiveMCMCHandler ( )

virtual

Destructor.

Definition at line 23 of file AdaptiveMCMCHandler.cpp.

                                           {
 // ********************************************
   if(adaptive_covariance != nullptr) {
     delete adaptive_covariance;
   }
 }

Member Function Documentation

◆ AdaptionUpdate()

bool AdaptiveMCMCHandler::AdaptionUpdate ( ) const

To be fair not a clue...

Definition at line 427 of file AdaptiveMCMCHandler.cpp.

                                                {
 // ********************************************
   if(total_steps <= start_adaptive_throw) return true;
   else return false;
 }

◆ CalculateRobbinsMonroStepLength()

void AdaptiveMCMCHandler::CalculateRobbinsMonroStepLength ( )

private

Calculate the constant step length for Robbins-Monro adaption.

Firstly we need to calculate the alpha value, this is in some sense "optimal"

Now we can calculate the scale factor

Definition at line 515 of file AdaptiveMCMCHandler.cpp.

                                                           {
 // ********************************************
   /* 
     Obtains the constant "step scale" for Robbins-Monro adaption
     for simplicity and because it has the degrees of freedom "baked in"
     it will be same across ALL blocks. 
   */
  
   double alpha = -ROOT::Math::normal_quantile(target_acceptance/2, 1.0); 
  
   int non_fixed_pars = GetNumParams()-GetNFixed();
  
   if(non_fixed_pars == 0){
     MACH3LOG_ERROR("Number of non_fixed_pars is 0, have you fixed all parameters");
     MACH3LOG_ERROR("This will cause division by 0 and crash ()");
     throw MaCh3Exception(__FILE__, __LINE__);
   }
   c_robbins_monro = (1- 1/non_fixed_pars)*std::sqrt(TMath::Pi()*2 * std::exp(alpha*alpha));
   c_robbins_monro += 1/(non_fixed_pars*target_acceptance*(1-target_acceptance));
  
   MACH3LOG_INFO("Robbins-Monro scale factor set to {}", c_robbins_monro);
 }

◆ CheckMatrixValidityForAdaption()

void AdaptiveMCMCHandler::CheckMatrixValidityForAdaption ( const TMatrixDSym * covMatrix ) const

Check if there are structures in matrix that could result in failure of adaption fits. this mostly cover cases of highly correlated block which seem to results in adaption failure.

Parameters

covMatrix covariance matrix which we convert into correlation matrix

Definition at line 449 of file AdaptiveMCMCHandler.cpp.

                                                                                            {
 // ********************************************
   int n = covMatrix->GetNrows();
   std::vector<std::vector<double>> corrMatrix(n, std::vector<double>(n));
  
   // Extract standard deviations from the diagonal
   std::vector<double> stdDevs(n);
   for (int i = 0; i < n; ++i) {
     stdDevs[i] = std::sqrt((*covMatrix)(i, i));
   }
  
   // Compute correlation for each element
   for (int i = 0; i < n; ++i) {
     for (int j = 0; j < n; ++j) {
       double cov = (*covMatrix)(i, j);
       double corr = cov / (stdDevs[i] * stdDevs[j]);
       corrMatrix[i][j] = corr;
     }
   }
  
   // KS: These numbers are completely arbitrary
   constexpr double NumberOfParametersThreshold = 5;
   constexpr double CorrelationThreshold = 0.5;
  
   bool FlagMessage = false;
   // For each parameter, count how many others it is correlated with
   for (int i = 0; i < n; ++i) {
     if (IsFixed(i)) {
       continue;
     }
     int block_i = adapt_block_matrix_indices[i];
     int correlatedCount = 0;
     std::vector<int> correlatedWith;
  
     for (int j = 0; j < n; ++j) {
       if (i == j || IsFixed(j) || adapt_block_matrix_indices[j] != block_i) {
         continue;
       }
       if (corrMatrix[i][j] > CorrelationThreshold) {
         correlatedCount++;
         correlatedWith.push_back(j);
       }
     }
  
     if (correlatedCount > NumberOfParametersThreshold) {
       MACH3LOG_WARN("Parameter {} ({}) is highly correlated with {} other parameters (above threshold {}).", _ParamNames->at(i), i, correlatedCount, CorrelationThreshold);
       MACH3LOG_WARN("Correlated with parameters: {}.", fmt::join(correlatedWith, ", "));
       FlagMessage = true;
     }
   }
   if(FlagMessage){
     MACH3LOG_WARN("Found highly correlated block, adaptive may not work as intended, proceed with caution");
     MACH3LOG_WARN("You can consider defying so called Adaption Block to not update this part of matrix");
   }
 }

◆ CreateNewAdaptiveCovariance()

void AdaptiveMCMCHandler::CreateNewAdaptiveCovariance ( )

If we don't have a covariance matrix to start from for adaptive tune we need to make one!

Definition at line 143 of file AdaptiveMCMCHandler.cpp.

                                                       {
 // ********************************************
   adaptive_covariance = new TMatrixDSym(GetNumParams());
   adaptive_covariance->Zero();
   par_means = std::vector<double>(GetNumParams(), 0);
 }

◆ CurrVal()

double AdaptiveMCMCHandler::CurrVal ( const int par_index ) const

Get Current value of parameter.

HW Implemented as its own method to allow for different behaviour in the future

Definition at line 507 of file AdaptiveMCMCHandler.cpp.

                                                              {
 // ********************************************
   return (*_fCurrVal)[par_index];
 }

◆ GetAdaptionScale()

double AdaptiveMCMCHandler::GetAdaptionScale ( ) const

inline

Get the current adaption scale.

Definition at line 216 of file AdaptiveMCMCHandler.h.

   {
     return adaption_scale;
   }

◆ GetAdaptiveCovariance()

const TMatrixDSym* AdaptiveMCMCHandler::GetAdaptiveCovariance ( ) const

inline

Return matrix used for adaption.

Definition at line 198 of file AdaptiveMCMCHandler.h.

   {
     return adaptive_covariance;
   }

◆ GetNFixed()

int AdaptiveMCMCHandler::GetNFixed ( ) const

inline

Get the number of fixed parameters.

Definition at line 138 of file AdaptiveMCMCHandler.h.

   {
     if (!_fFixedPars)
     {
       return 0;
     }
     int n_fixed = 0;
     for (int i = 0; i < static_cast<int>(_fFixedPars->size()); i++)
     {
       if ((*_fFixedPars)[i] < 0)
       {
         n_fixed++;
       }
     }
     return n_fixed;
   }

◆ GetNumParams()

int AdaptiveMCMCHandler::GetNumParams ( ) const

inline

Get the current values of the parameters.

Definition at line 156 of file AdaptiveMCMCHandler.h.

   {
     return static_cast<int>(_fCurrVal->size());
   }

◆ GetOutFileName()

std::string AdaptiveMCMCHandler::GetOutFileName ( ) const

inline

Get Name of Output File.

Definition at line 210 of file AdaptiveMCMCHandler.h.

   {
     return output_file_name;
   }

◆ GetParameterMeans()

std::vector<double> AdaptiveMCMCHandler::GetParameterMeans ( ) const

inline

Get the parameter means used in the adaptive handler.

Definition at line 204 of file AdaptiveMCMCHandler.h.

   {
     return par_means;
   }

◆ GetParIndex()

int AdaptiveMCMCHandler::GetParIndex ( const std::string & name ) const

inline

Get the index of the parameter given its name.

Parameters

name	Name of the parameter

Returns: Index of the parameter

Definition at line 96 of file AdaptiveMCMCHandler.h.

                                                {
     // HH: Adapted from ParameterHandlerBase.cpp
     int index = M3::_BAD_INT_;
     for (size_t i = 0; i < _ParamNames->size(); i++) {
       if(name == _ParamNames->at(i)) {
         index = static_cast<int>(i);
         break;
       }
     }
     return index;
   }

◆ GetParIndicesFromNames()

std::vector<int> AdaptiveMCMCHandler::GetParIndicesFromNames ( const std::vector< std::string > & param_names ) const

inline

Convert vector of parameter names to indices.

Parameters

param_names Vector of parameter names

Returns: Vector of parameter indices

Definition at line 111 of file AdaptiveMCMCHandler.h.

                                                                                        {
     std::vector<int> indices;
     for (const auto& name : param_names) {
       int index = GetParIndex(name);
       if (index != M3::_BAD_INT_) {
         indices.push_back(index);
       } else {
         MACH3LOG_ERROR("Parameter name {} not found in parameter list", name);
         throw MaCh3Exception(__FILE__, __LINE__);
       }
     }
     return indices;
   }

◆ GetTotalSteps()

int AdaptiveMCMCHandler::GetTotalSteps ( ) const

inline

Get Total Number of Steps.

Definition at line 175 of file AdaptiveMCMCHandler.h.

   {
     return total_steps;
   }

◆ GetUseRobbinsMonro()

bool AdaptiveMCMCHandler::GetUseRobbinsMonro ( ) const

inline

Use Robbins-Monro approach?

Definition at line 222 of file AdaptiveMCMCHandler.h.

   {
     return use_robbins_monro;
   }

◆ IncrementAcceptedSteps()

void AdaptiveMCMCHandler::IncrementAcceptedSteps ( )

inline

Definition at line 192 of file AdaptiveMCMCHandler.h.

   {
     prev_step_accepted = true;
   }

◆ IncrementNSteps()

void AdaptiveMCMCHandler::IncrementNSteps ( )

inline

Increase by one number of total steps.

Definition at line 187 of file AdaptiveMCMCHandler.h.

   {
     total_steps++;
   }

◆ IndivStepScaleAdapt()

bool AdaptiveMCMCHandler::IndivStepScaleAdapt ( ) const

Tell whether we want reset step scale or not.

Definition at line 400 of file AdaptiveMCMCHandler.cpp.

                                                     {
 // ********************************************
   if(total_steps == start_adaptive_throw) return true;
   else return false;
 }

◆ InitFromConfig()

bool AdaptiveMCMCHandler::InitFromConfig	(	const YAML::Node &	adapt_manager,
		const std::string &	matrix_name_str,
		const std::vector< std::string > *	parameter_names,
		const std::vector< double > *	parameters,
		const std::vector< double > *	fixed,
		const std::vector< bool > *	param_skip_adapt_flags,
		const TMatrixDSym *	throwMatrix,
		const double	initial_scale_
	)

Read initial values from config file.

Parameters

adapt_manager	YAML node containing the configuration (`AdaptionOptions`).
matrix_name_str	Name of the covariance matrix block to configure.
parameters	Pointer to a vector of parameter values (nominal values).
fixed	Pointer to a vector of fixed parameter values.

Returns: True if adaptive MCMC configuration was successfully initialized, false otherwise.

HW: This is technically wrong, should be across all systematics but will be addressed in a later PR

Definition at line 31 of file AdaptiveMCMCHandler.cpp.

                                                                      {
 // ********************************************
   /*
    * HW: Idea is that adaption can simply read the YAML config
    * Options :
    *         External Info:
    * UseExternalMatrix [bool]     :    Use an external matrix
    * ExternalMatrixFileName [str] :    Name of file containing external info
    * ExternalMatrixName [str]     :    Name of external Matrix
    * ExternalMeansName [str]      :    Name of external means vector [for updates]
    *
    *         General Info:
    * DoAdaption [bool]            :    Do we want to do adaption?
    * AdaptionStartThrow [int]     :    Step we start throwing adaptive matrix from
    * AdaptionEndUpdate [int]      :    Step we stop updating adaptive matrix
    * AdaptionStartUpdate [int]    :    Do we skip the first N steps?
    * AdaptionUpdateStep [int]     :    Number of steps between matrix updates
    * AdaptionSaveNIterations [int]:    You don't have to save every adaptive stage so decide how often you want to save
    * Adaption blocks [vector<vector<int>>] : Splits the throw matrix into several block matrices
    * OuputFileName [std::string]  :    Name of the file that the adaptive matrices will be saved into
    */
  
   // setAdaptionDefaults();
   if(!adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]) {
     MACH3LOG_WARN("Adaptive Settings not found for {}, this is fine if you don't want adaptive MCMC", matrix_name_str);
     return false;
   }
  
   // We"re going to grab this info from the YAML manager
   if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["DoAdaption"], false, __FILE__ , __LINE__)) {
     MACH3LOG_WARN("Not using adaption for {}", matrix_name_str);
     return false;
   }
  
   if(!CheckNodeExists(adapt_manager, "AdaptionOptions", "Settings", "OutputFileName")) {
     MACH3LOG_ERROR("No OutputFileName specified in AdaptionOptions::Settings into your config file");
     MACH3LOG_ERROR("This is required if you are using adaptive MCMC");
     throw MaCh3Exception(__FILE__, __LINE__);
   }
   auto AdaptSettings = adapt_manager["AdaptionOptions"]["Settings"];
   start_adaptive_throw  = GetFromManager<int>(AdaptSettings["StartThrow"], 10, __FILE__ , __LINE__);
   start_adaptive_update = GetFromManager<int>(AdaptSettings["StartUpdate"], 0, __FILE__ , __LINE__);
   end_adaptive_update   = GetFromManager<int>(AdaptSettings["EndUpdate"], 10000, __FILE__ , __LINE__);
   adaptive_update_step  = GetFromManager<int>(AdaptSettings["UpdateStep"], 100, __FILE__ , __LINE__);
   adaptive_save_n_iterations  = GetFromManager<int>(AdaptSettings["SaveNIterations"], -1, __FILE__ , __LINE__);
   output_file_name = GetFromManager<std::string>(AdaptSettings["OutputFileName"], "", __FILE__ , __LINE__);
  
   // Check for Robbins-Monro adaption
   use_robbins_monro = GetFromManager<bool>(AdaptSettings["UseRobbinsMonro"], false, __FILE__ , __LINE__);
   target_acceptance = GetFromManager<double>(AdaptSettings["TargetAcceptance"], 0.234, __FILE__ , __LINE__);
   if (target_acceptance <= 0 || target_acceptance >= 1) {
     MACH3LOG_ERROR("Target acceptance must be in (0,1), got {}", target_acceptance);
     throw MaCh3Exception(__FILE__, __LINE__);
   }
  
   acceptance_rate_batch_size = GetFromManager<int>(AdaptSettings["AcceptanceRateBatchSize"], 10000, __FILE__ , __LINE__);
   total_rm_restarts = GetFromManager<int>(AdaptSettings["TotalRobbinsMonroRestarts"], 0, __FILE__ , __LINE__);
   n_rm_restarts = 0;
  
   prev_step_accepted = false; 
  
   // We also want to check for "blocks" by default all parameters "know" about each other
   // but we can split the matrix into independent block matrices
   SetParams(parameters);
   SetFixed(fixed);
  
   _ParamNames = parameter_names;
  
   initial_throw_matrix = static_cast<TMatrixDSym*>(throwMatrix->Clone());
  
   initial_scale = initial_scale_;
  
   if(use_robbins_monro) {
     MACH3LOG_INFO("Using Robbins-Monro for adaptive step size tuning with target acceptance {}", target_acceptance);
     // Now we need some thing
     CalculateRobbinsMonroStepLength();
   }
  
   // We'll use the same start scale for Robbins-Monro as standard adaption
   adaption_scale = 2.38/std::sqrt(GetNumParams()); 
  
   auto AdaptOptions = adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str];
   // HH: added ability to define blocks by parameter names
   bool matrix_blocks_by_name = GetFromManager<bool>(AdaptOptions["MatrixBlocksByName"], false, __FILE__ , __LINE__);
   std::vector<std::vector<int>> matrix_blocks;
   // HH: read blocks by names and convert to indices
   if (matrix_blocks_by_name) {
     auto matrix_blocks_input = GetFromManager<std::vector<std::vector<std::string>>>(AdaptOptions["MatrixBlocks"], {{}}, __FILE__ , __LINE__);
     // Now we need to convert to indices
     for (const auto& block : matrix_blocks_input) {
       auto block_indices = GetParIndicesFromNames(block);
       matrix_blocks.push_back(block_indices);
     }
   } else {
     matrix_blocks = GetFromManager<std::vector<std::vector<int>>>(AdaptOptions["MatrixBlocks"], {{}}, __FILE__ , __LINE__);
   }
   SetAdaptiveBlocks(matrix_blocks);
  
   // set parameters to skip during adaption
   _param_skip_adapt_flags = param_skip_adapt_flags;
  
   return true;
 }

◆ IsFixed()

bool AdaptiveMCMCHandler::IsFixed ( const int ipar ) const

inline

Check if a parameter is fixed.

Definition at line 162 of file AdaptiveMCMCHandler.h.

   {
     if (!_fFixedPars)
     {
       return false;
     }
     return ((*_fFixedPars)[ipar] < 0);
   }

◆ Print()

void AdaptiveMCMCHandler::Print ( ) const

Print all class members.

Definition at line 434 of file AdaptiveMCMCHandler.cpp.

                                       {
 // ********************************************
   MACH3LOG_INFO("Adaptive MCMC Info:");
   MACH3LOG_INFO("Throwing from New Matrix from Step : {}", start_adaptive_throw);
   MACH3LOG_INFO("Adaption Matrix Start Update       : {}", start_adaptive_update);
   MACH3LOG_INFO("Adaption Matrix Ending Updates     : {}", end_adaptive_update);
   MACH3LOG_INFO("Steps Between Updates              : {}", adaptive_update_step);
   MACH3LOG_INFO("Saving matrices to file            : {}", output_file_name);
   MACH3LOG_INFO("Will only save every {} iterations"     , adaptive_save_n_iterations); 
   if(use_robbins_monro) {
     MACH3LOG_INFO("Using Robbins-Monro for adaptive step size tuning with target acceptance {}", target_acceptance);
   }
 }

◆ SaveAdaptiveToFile()

void AdaptiveMCMCHandler::SaveAdaptiveToFile	(	const std::string &	outFileName,
		const std::string &	systematicName,
		const bool	is_final = `false`
	)

HW: Save adaptive throw matrix to file.

Definition at line 193 of file AdaptiveMCMCHandler.cpp.

                                                                                              {
 // ********************************************
   // Skip saving if adaptive_save_n_iterations is negative,
   // unless this is the final iteration (is_final overrides the condition)
   if (adaptive_save_n_iterations < 0 && !is_final) return;
   if (is_final ||
     (adaptive_save_n_iterations > 0 &&
     ((total_steps - start_adaptive_throw) / adaptive_update_step) % adaptive_save_n_iterations == 0)) {
     TFile *outFile = new TFile(outFileName.c_str(), "UPDATE");
     if (outFile->IsZombie()) {
       MACH3LOG_ERROR("Couldn't find {}", outFileName);
       throw MaCh3Exception(__FILE__, __LINE__);
     }
  
     TVectorD *outMeanVec = new TVectorD(int(par_means.size()));
     for (int i = 0; i < int(par_means.size()); i++) {
       (*outMeanVec)(i) = par_means[i];
     }
  
     std::string adaptive_cov_name = systematicName + "_posfit_matrix";
     std::string mean_vec_name = systematicName + "_mean_vec";
     if (!is_final) {
       // Some string to make the name of the saved adaptive matrix clear
       std::string total_steps_str =
           std::to_string(total_steps);
       std::string syst_name_str = systematicName;
       adaptive_cov_name =
           total_steps_str + '_' + syst_name_str + std::string("_throw_matrix");
       mean_vec_name =
           total_steps_str + '_' + syst_name_str + std::string("_mean_vec");
     }
  
     outFile->cd();
     auto adaptive_to_save = static_cast<TMatrixDSym*>(adaptive_covariance->Clone());
     // Multiply by scale
     // HH: Only scale parameters not being skipped
     for (int i = 0; i < GetNumParams(); ++i) {
       for (int j = 0; j <= i; ++j) {
         // If both parameters are skipped, scale by initial scale
         if ((*_param_skip_adapt_flags)[i] && (*_param_skip_adapt_flags)[j]) {
           (*adaptive_to_save)(i, j) *= initial_scale * initial_scale;
           if (i != j) {
             (*adaptive_to_save)(j, i) *= initial_scale * initial_scale;
           }
         }
         else if (!(*_param_skip_adapt_flags)[i] && !(*_param_skip_adapt_flags)[j]) {
           (*adaptive_to_save)(i, j) *= GetAdaptionScale() * GetAdaptionScale();
           if (i != j) {
             (*adaptive_to_save)(j, i) *= GetAdaptionScale() * GetAdaptionScale();
           }
         }
         else {
           // Not too sure what to do here, as the correlation should be zero 
           // Scale by both factors as a compromise
           (*adaptive_to_save)(i, j) *= GetAdaptionScale() * initial_scale;
           if (i != j) {
             (*adaptive_to_save)(j, i) *= GetAdaptionScale() * initial_scale;
           }
         }
       }
     }
  
     adaptive_to_save->Write(adaptive_cov_name.c_str());
     outMeanVec->Write(mean_vec_name.c_str());
  
     // HH: Also save the initial throw matrix for reference
     if (!initial_throw_matrix_saved) {
       std::string initial_cov_name = "0_" + systematicName + std::string("_initial_throw_matrix");
       auto initial_throw_matrix_to_save = static_cast<TMatrixDSym*>(initial_throw_matrix->Clone());
       *(initial_throw_matrix_to_save) *= initial_scale * initial_scale;
       initial_throw_matrix_saved = true;
       initial_throw_matrix_to_save->Write(initial_cov_name.c_str());
       delete initial_throw_matrix_to_save;
     }
  
     outFile->Close();
     delete adaptive_to_save;
     delete outMeanVec;
     delete outFile;
   }
 }

◆ SetAdaptiveBlocks()

void AdaptiveMCMCHandler::SetAdaptiveBlocks ( const std::vector< std::vector< int >> & block_indices )

HW: sets adaptive block matrix.

Parameters

block_indices Values for sub-matrix blocks

Definition at line 151 of file AdaptiveMCMCHandler.cpp.

                                                                                           {
 // ********************************************
   /*
    *   In order to adapt efficient we want to setup our throw matrix to be a serious of block-diagonal (ish) matrices
    *
    *   To do this we set sub-block in the config by parameter index. For example having
    *   [[0,4],[4, 6]] in your config will set up two blocks one with all indices 0<=i<4 and the other with 4<=i<6
    */
   // Set up block regions
   adapt_block_matrix_indices = std::vector<int>(GetNumParams(), 0);
  
   // Should also make a matrix of block sizes
   adapt_block_sizes = std::vector<int>(block_indices.size()+1, 0);
   adapt_block_sizes[0] = GetNumParams();
  
   if(block_indices.size()==0 || block_indices[0].size()==0) return;
  
   int block_size = static_cast<int>(block_indices.size());
   // Now we loop over our blocks
   for(int iblock=0; iblock < block_size; iblock++){
     // Loop over blocks in the block
     int sub_block_size = static_cast<int>(block_indices.size()-1);
     for(int isubblock=0; isubblock < sub_block_size ; isubblock+=2){
       int block_lb = block_indices[iblock][isubblock];
       int block_ub = block_indices[iblock][isubblock+1];
  
       if(block_lb > GetNumParams() || block_ub > GetNumParams()){
         MACH3LOG_ERROR("Cannot set matrix block with edges {}, {} for matrix of size {}",
                        block_lb, block_ub, GetNumParams());
         throw MaCh3Exception(__FILE__, __LINE__);;
       }
       for(int ipar = block_lb; ipar < block_ub; ipar++){
         adapt_block_matrix_indices[ipar] = iblock+1;
         adapt_block_sizes[iblock+1] += 1;
         adapt_block_sizes[0] -= 1;
       }
     }
   }
 }

◆ SetFixed()

void AdaptiveMCMCHandler::SetFixed ( const std::vector< double > * fix )

inline

Set the fixed parameters.

Definition at line 132 of file AdaptiveMCMCHandler.h.

   {
     _fFixedPars = fix;
   }

◆ SetParams()

void AdaptiveMCMCHandler::SetParams ( const std::vector< double > * params )

inline

Set the current values of the parameters.

Definition at line 126 of file AdaptiveMCMCHandler.h.

   {
     _fCurrVal = params;
   }

◆ SetThrowMatrixFromFile()

void AdaptiveMCMCHandler::SetThrowMatrixFromFile	(	const std::string &	matrix_file_name,
		const std::string &	matrix_name,
		const std::string &	means_name,
		bool &	use_adaptive
	)

sets throw matrix from a file

Parameters

matrix_file_name	name of file matrix lives in
matrix_name	name of matrix in file
means_name	name of means vec in file

Definition at line 279 of file AdaptiveMCMCHandler.cpp.

                                                                      {
 // ********************************************
   // Lets you set the throw matrix externally
   // Open file
   auto matrix_file = std::make_unique<TFile>(matrix_file_name.c_str());
   use_adaptive = true;
  
   if(matrix_file->IsZombie()){
     MACH3LOG_ERROR("Couldn't find {}", matrix_file_name);
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
  
   // Next we grab our matrix
   adaptive_covariance = static_cast<TMatrixDSym*>(matrix_file->Get(matrix_name.c_str()));
   if(!adaptive_covariance){
     MACH3LOG_ERROR("Couldn't find {} in {}", matrix_name, matrix_file_name);
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
  
   // Finally we grab the means vector
   TVectorD* means_vector = static_cast<TVectorD*>(matrix_file->Get(means_name.c_str()));
  
   // This is fine to not exist!
   if(means_vector){
     // Yay our vector exists! Let's loop and fill it
     // Should check this is done
     if(means_vector->GetNrows() != GetNumParams()){
       MACH3LOG_ERROR("External means vec size ({}) != matrix size ({})", means_vector->GetNrows(), GetNumParams());
       throw MaCh3Exception(__FILE__, __LINE__);
     }
  
     par_means = std::vector<double>(GetNumParams());
     for(int i = 0; i < GetNumParams(); i++){
       par_means[i] = (*means_vector)(i);
     }
     MACH3LOG_INFO("Found Means in External File, Will be able to adapt");
   }
   // Totally fine if it doesn't exist, we just can't do adaption
   else{
     // We don't need a means vector, set the adaption=false
     MACH3LOG_WARN("Cannot find means vector in {}, therefore I will not be able to adapt!", matrix_file_name);
     use_adaptive = false;
   }
  
   matrix_file->Close();
   MACH3LOG_INFO("Set up matrix from external file");
 }

◆ SetTotalSteps()

void AdaptiveMCMCHandler::SetTotalSteps ( const int nsteps )

inline

Change Total Number of Steps to new value.

Definition at line 181 of file AdaptiveMCMCHandler.h.

   {
     total_steps = nsteps;
   }

◆ SkipAdaption()

bool AdaptiveMCMCHandler::SkipAdaption ( ) const

Tell if we are Skipping Adaption.

Definition at line 419 of file AdaptiveMCMCHandler.cpp.

                                              {
 // ********************************************
   if(total_steps > end_adaptive_update ||
     total_steps< start_adaptive_update) return true;
   else return false;
 }

◆ UpdateAdaptiveCovariance()

void AdaptiveMCMCHandler::UpdateAdaptiveCovariance ( )

Method to update adaptive MCMC [15].

Parameters

_fCurrVal Value of each parameter necessary for updating throw matrix

Definition at line 331 of file AdaptiveMCMCHandler.cpp.

                                                    {
 // ********************************************
   std::vector<double> par_means_prev = par_means;
   int steps_post_burn = total_steps - start_adaptive_update;
  
   prev_step_accepted = false;
  
   // Step 1: Update means and compute deviations
   for (int i = 0; i < GetNumParams(); ++i) {
     if (IsFixed(i)) continue;
  
     par_means[i] =  (CurrVal(i) + par_means_prev[i]*steps_post_burn)/(steps_post_burn+1);
     // Left over from cyclic means
   }
  
   // Step 2: Update covariance
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif
   for (int i = 0; i < GetNumParams(); ++i) {
     if (IsFixed(i)) {
       (*adaptive_covariance)(i, i) = 1.0;
       continue;
     }
  
     int block_i = adapt_block_matrix_indices[i];
  
     for (int j = 0; j <= i; ++j) {
       if (IsFixed(j) || adapt_block_matrix_indices[j] != block_i) {
         (*adaptive_covariance)(i, j) = 0.0;
         (*adaptive_covariance)(j, i) = 0.0;
         continue;
       }
  
       double cov_prev = (*adaptive_covariance)(i, j);
       double cov_updated = 0.0;
  
       // HH: Check if either parameter is skipped
       // If parameter i and j are both skipped
       if ((*_param_skip_adapt_flags)[i] && (*_param_skip_adapt_flags)[j]) {
         // Set covariance to initial throw matrix value
         (*adaptive_covariance)(i, j) = (*initial_throw_matrix)(i, j);
         (*adaptive_covariance)(j, i) = (*initial_throw_matrix)(j, i);
         continue;
       }
       // If only one parameter is skipped we need to make sure correlation is zero 
       // because we don't want to change one parameter while keeping the other fixed 
       // as this would lead to penalty terms in the prior
       else if ((*_param_skip_adapt_flags)[i] || (*_param_skip_adapt_flags)[j]) {
         continue;
       } // otherwise adapt as normal
  
       if (steps_post_burn > 0) {
       // Haario-style update
         double cov_t = cov_prev * (steps_post_burn - 1) / steps_post_burn;
         double prev_means_t = steps_post_burn * par_means_prev[i] * par_means_prev[j];
         double curr_means_t = (steps_post_burn + 1) * par_means[i] * par_means[j];
         double curr_step_t = CurrVal(i) * CurrVal(j);
  
         cov_updated = cov_t + (prev_means_t - curr_means_t + curr_step_t) / steps_post_burn;
       }
  
       (*adaptive_covariance)(i, j) = cov_updated;
       (*adaptive_covariance)(j, i) = cov_updated;
     }
   }
 }

◆ UpdateMatrixAdapt()

bool AdaptiveMCMCHandler::UpdateMatrixAdapt ( )

Tell whether matrix should be updated.

Definition at line 407 of file AdaptiveMCMCHandler.cpp.

                                             {
 // ********************************************
   if(total_steps >= start_adaptive_throw &&
     // Check whether the number of steps is divisible by the adaptive update step
     // e.g. if adaptive_update_step = 1000 and (total_step - start_adpative_throw) is 5000 then this is true
     (total_steps - start_adaptive_throw)%adaptive_update_step == 0) {
     return true;
   } 
   else return false;
 }

◆ UpdateRobbinsMonroScale()

void AdaptiveMCMCHandler::UpdateRobbinsMonroScale ( )

Update the scale factor for Robbins-Monro adaption.

This allows to move towards a STABLE number of steps in the batch

Update the scale factor for Robbins-Monro adaption [200 is arbitrary for early stability but motivated by paper]

Now we either increase or decrease the scale

Definition at line 541 of file AdaptiveMCMCHandler.cpp.

                                                  {
 // ********************************************
   /* 
   Update the scale factor using Robbins-Monro. 
   TLDR: If acceptance rate is too high, scale factor goes up, if too low goes down
         will pretty rapidly converge to the right value.
   */
   if ((total_steps>0) && (total_steps % acceptance_rate_batch_size == 0)) {
     n_rm_restarts++;
   }
  
   int batch_step;
   if (n_rm_restarts > total_rm_restarts) {
     batch_step = total_steps;
   }
   else{
     batch_step = total_steps % acceptance_rate_batch_size;
   }
  
   int non_fixed_pars = GetNumParams()-GetNFixed();
  
   // double root_scale = std::sqrt(adaption_scale);
  
   double scale_factor = adaption_scale * c_robbins_monro / std::max(200.0, static_cast<double>(batch_step) / non_fixed_pars);
  
   if(prev_step_accepted){
     adaption_scale += (1 - target_acceptance) * scale_factor;
   }
   else{
     adaption_scale -= target_acceptance * scale_factor;
   }
 }

Member Data Documentation

◆ _fCurrVal

const std::vector<double>* AdaptiveMCMCHandler::_fCurrVal

private

Current values of parameters.

Definition at line 280 of file AdaptiveMCMCHandler.h.

◆ _fFixedPars

const std::vector<double>* AdaptiveMCMCHandler::_fFixedPars

private

Vector of fixed parameters.

Definition at line 277 of file AdaptiveMCMCHandler.h.

◆ _param_skip_adapt_flags

const std::vector<bool>* AdaptiveMCMCHandler::_param_skip_adapt_flags

private

Parameters to skip during adaption.

Definition at line 307 of file AdaptiveMCMCHandler.h.

◆ _ParamNames

const std::vector<std::string>* AdaptiveMCMCHandler::_ParamNames

private

Parameter names.

Definition at line 283 of file AdaptiveMCMCHandler.h.

◆ acceptance_rate_batch_size

int AdaptiveMCMCHandler::acceptance_rate_batch_size

private

Acceptance rate in the current batch.

Definition at line 286 of file AdaptiveMCMCHandler.h.

◆ adapt_block_matrix_indices

std::vector<int> AdaptiveMCMCHandler::adapt_block_matrix_indices

private

Indices for block-matrix adaption.

Definition at line 255 of file AdaptiveMCMCHandler.h.

◆ adapt_block_sizes

std::vector<int> AdaptiveMCMCHandler::adapt_block_sizes

private

Size of blocks for adaption.

Definition at line 258 of file AdaptiveMCMCHandler.h.

◆ adaption_scale

double AdaptiveMCMCHandler::adaption_scale

private

Scaling factor.

Definition at line 271 of file AdaptiveMCMCHandler.h.

◆ adaptive_covariance

TMatrixDSym* AdaptiveMCMCHandler::adaptive_covariance

private

Full adaptive covariance matrix.

Definition at line 265 of file AdaptiveMCMCHandler.h.

◆ adaptive_save_n_iterations

int AdaptiveMCMCHandler::adaptive_save_n_iterations

private

If you don't want to save every adaption then you can specify this here

Definition at line 249 of file AdaptiveMCMCHandler.h.

◆ adaptive_update_step

int AdaptiveMCMCHandler::adaptive_update_step

private

Steps between changing throw matrix.

Definition at line 245 of file AdaptiveMCMCHandler.h.

◆ c_robbins_monro

double AdaptiveMCMCHandler::c_robbins_monro

private

Constant "step scaling" factor for Robbins-Monro.

Definition at line 295 of file AdaptiveMCMCHandler.h.

◆ end_adaptive_update

int AdaptiveMCMCHandler::end_adaptive_update

private

Steps between changing throw matrix.

Definition at line 242 of file AdaptiveMCMCHandler.h.

◆ initial_scale

double AdaptiveMCMCHandler::initial_scale

private

Initial scaling factor.

Definition at line 274 of file AdaptiveMCMCHandler.h.

◆ initial_throw_matrix

const TMatrixDSym* AdaptiveMCMCHandler::initial_throw_matrix

private

Initial throw matrix.

Definition at line 310 of file AdaptiveMCMCHandler.h.

◆ initial_throw_matrix_saved

bool AdaptiveMCMCHandler::initial_throw_matrix_saved = false

private

Flag to check if initial throw matrix has been saved.

Definition at line 313 of file AdaptiveMCMCHandler.h.

◆ n_rm_restarts

int AdaptiveMCMCHandler::n_rm_restarts

private

Number of restarts for Robbins Monro (so far)

Definition at line 298 of file AdaptiveMCMCHandler.h.

◆ output_file_name

std::string AdaptiveMCMCHandler::output_file_name

private

Name of the file to save the adaptive matrices into.

Definition at line 252 of file AdaptiveMCMCHandler.h.

◆ par_means

std::vector<double> AdaptiveMCMCHandler::par_means

private

Mean values for all parameters.

Definition at line 262 of file AdaptiveMCMCHandler.h.

◆ prev_step_accepted

bool AdaptiveMCMCHandler::prev_step_accepted

private

Need to keep track of whether previous step was accepted for RM.

Definition at line 304 of file AdaptiveMCMCHandler.h.

◆ start_adaptive_throw

int AdaptiveMCMCHandler::start_adaptive_throw

private

Meta variables related to adaption run time When do we start throwing

Definition at line 236 of file AdaptiveMCMCHandler.h.

◆ start_adaptive_update

int AdaptiveMCMCHandler::start_adaptive_update

private

When do we stop update the adaptive matrix.

Definition at line 239 of file AdaptiveMCMCHandler.h.

◆ target_acceptance

double AdaptiveMCMCHandler::target_acceptance

private

Target acceptance rate for Robbins Monro.

Definition at line 292 of file AdaptiveMCMCHandler.h.

◆ total_rm_restarts

int AdaptiveMCMCHandler::total_rm_restarts

private

Total number of restarts ALLOWED for Robbins Monro.

Definition at line 301 of file AdaptiveMCMCHandler.h.

◆ total_steps

int AdaptiveMCMCHandler::total_steps

private

Total number of MCMC steps.

Definition at line 268 of file AdaptiveMCMCHandler.h.

◆ use_robbins_monro

bool AdaptiveMCMCHandler::use_robbins_monro

private

Use Robbins Monro https://arxiv.org/pdf/1006.3690.

Definition at line 289 of file AdaptiveMCMCHandler.h.

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

Parameters/AdaptiveMCMCHandler.h
Parameters/AdaptiveMCMCHandler.cpp

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ AdaptiveMCMCHandler()

◆ ~AdaptiveMCMCHandler()

Member Function Documentation

◆ AdaptionUpdate()

◆ CalculateRobbinsMonroStepLength()

◆ CheckMatrixValidityForAdaption()

◆ CreateNewAdaptiveCovariance()

◆ CurrVal()

◆ GetAdaptionScale()

◆ GetAdaptiveCovariance()

◆ GetNFixed()

◆ GetNumParams()

◆ GetOutFileName()

◆ GetParameterMeans()

◆ GetParIndex()

◆ GetParIndicesFromNames()

◆ GetTotalSteps()

◆ GetUseRobbinsMonro()

◆ IncrementAcceptedSteps()

◆ IncrementNSteps()

◆ IndivStepScaleAdapt()

◆ InitFromConfig()

◆ IsFixed()

◆ Print()

◆ SaveAdaptiveToFile()

◆ SetAdaptiveBlocks()

◆ SetFixed()

◆ SetParams()

◆ SetThrowMatrixFromFile()

◆ SetTotalSteps()

◆ SkipAdaption()

◆ UpdateAdaptiveCovariance()

◆ UpdateMatrixAdapt()

◆ UpdateRobbinsMonroScale()

Member Data Documentation

◆ _fCurrVal

◆ _fFixedPars

◆ _param_skip_adapt_flags

◆ _ParamNames

◆ acceptance_rate_batch_size

◆ adapt_block_matrix_indices

◆ adapt_block_sizes

◆ adaption_scale

◆ adaptive_covariance

◆ adaptive_save_n_iterations

◆ adaptive_update_step

◆ c_robbins_monro

◆ end_adaptive_update

◆ initial_scale

◆ initial_throw_matrix

◆ initial_throw_matrix_saved

◆ n_rm_restarts

◆ output_file_name

◆ par_means

◆ prev_step_accepted

◆ start_adaptive_throw

◆ start_adaptive_update

◆ target_acceptance

◆ total_rm_restarts

◆ total_steps

◆ use_robbins_monro