ParameterHandlerBase Class Reference

Base class for handling systematic uncertainty parameters. More...

#include <Parameters/ParameterHandlerBase.h>

Inheritance diagram for ParameterHandlerBase:

Collaboration diagram for ParameterHandlerBase:

Public Member Functions
	ParameterHandlerBase (std::string name, std::string file, double threshold=-1, int FirstPCAdpar=-999, int LastPCAdpar=-999)
	"Usual" constructors from root file More...
	ParameterHandlerBase ()=default
virtual	~ParameterHandlerBase ()
	Destructor. More...
void	SetCovMatrix (TMatrixDSym *cov)
	Set covariance matrix. More...
void	SetName (const std::string &name)
	Set matrix name. More...
void	SetParName (const int i, const std::string &name)
	change parameter name More...
void	SetSingleParameter (const int parNo, const double parVal)
	Set value of single param to a given value. More...
void	SetPar (const int i, const double val)
	Set all the covariance matrix parameters to a user-defined value. More...
void	SetParCurrProp (const int i, const double val)
	Set current parameter value. More...
void	SetParProp (const int i, const double val)
	Set proposed parameter value. More...
void	SetParameters (const std::vector< double > &pars={})
	Set parameter values using vector, it has to have same size as covariance class. More...
void	SetFlatPrior (const int i, const bool eL)
	Set if parameter should have flat prior or not. More...
void	SetRandomThrow (const int i, const double rand)
	Set random value useful for debugging/CI. More...
double	GetRandomThrow (const int i) const
	Get random value useful for debugging/CI. More...
void	SetBranches (TTree &tree, const bool SaveProposal=false)
	set branches for output file More...
void	SetStepScale (const double scale, const bool verbose=true)
	Set global step scale for covariance object. More...
void	SetIndivStepScale (const int ParameterIndex, const double StepScale)
	DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constructors) More...
void	SetIndivStepScale (const std::vector< double > &stepscale)
	DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constructors) More...
void	SetPrintLength (const unsigned int PriLen)
	KS: In case someone really want to change this. More...
void	SaveUpdatedMatrixConfig ()
	KS: After step scale, prefit etc. value were modified save this modified config. More...
void	ThrowParameters ()
	Throw the parameters according to the covariance matrix. This shouldn't be used in MCMC code ase it can break Detailed Balance;. More...
void	RandomConfiguration ()
	Randomly throw the parameters in their 1 sigma range. More...
int	CheckBounds () const _noexcept_
	Check if parameters were proposed outside physical boundary. More...
double	CalcLikelihood () const _noexcept_
	Calc penalty term based on inverted covariance matrix. More...
virtual double	GetLikelihood ()
	Return CalcLikelihood if some params were thrown out of boundary return LARGE_LOGL More...
TMatrixDSym *	GetCovMatrix () const
	Return covariance matrix. More...
TMatrixDSym *	GetInvCovMatrix () const
	Return inverted covariance matrix. More...
double	GetInvCovMatrix (const int i, const int j) const
	Return inverted covariance matrix. More...
double	GetCorrThrows (const int i) const
	Return correlated throws. More...
bool	GetFlatPrior (const int i) const
	Get if param has flat prior or not. More...
std::string	GetName () const
	Get name of covariance. More...
std::string	GetParName (const int i) const
	Get name of parameter. More...
int	GetParIndex (const std::string &name) const
	Get index based on name. More...
std::string	GetParFancyName (const int i) const
	Get fancy name of the Parameter. More...
double	GetDiagonalError (const int i) const
	Get diagonal error for ith parameter. More...
double	GetError (const int i) const
	Get the error for the ith parameter. More...
void	ResetIndivStepScale ()
	Adaptive Step Tuning Stuff. More...
void	SetIndivStepScaleForSkippedAdaptParams ()
	Set individual step scale for parameters which are skipped during adaption to initial values. More...
void	InitialiseAdaption (const YAML::Node &adapt_manager)
	Initialise adaptive MCMC. More...
void	SaveAdaptiveToFile (const std::string &outFileName, const std::string &systematicName)
	Save adaptive throw matrix to file. More...
bool	GetDoAdaption () const
	Do we adapt or not. More...
void	SetThrowMatrix (const TMatrixDSym *cov)
	Use new throw matrix, used in adaptive MCMC. More...
void	SetSubThrowMatrix (int first_index, int last_index, TMatrixDSym const &subcov)
void	UpdateThrowMatrix (TMatrixDSym *cov)
	Replaces old throw matrix with new one. More...
void	SetNumberOfSteps (const int nsteps)
	Set number of MCMC step, when running adaptive MCMC it is updated with given frequency. We need number of steps to determine frequency. More...
TMatrixDSym *	GetThrowMatrix () const
	Get matrix used for step proposal. More...
double	GetThrowMatrix (const int i, const int j) const
	Get matrix used for step proposal. More...
TH2D *	GetCorrelationMatrix () const
	KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting. More...
const M3::float_t *	RetPointer (const int iParam) const
	DB Pointer return to param position. More...
const std::vector< M3::float_t > &	GetParPropVec ()
	Get a reference to the proposed parameter values Can be useful if you want to track these without having to copy values using getProposed() More...
int	GetNumParams () const
	Get total number of parameters. More...
std::vector< double >	GetPreFitValues () const
	Get the pre-fit values of the parameters. More...
std::vector< double >	GetProposed () const
	Get vector of all proposed parameter values. More...
M3::float_t	GetParProp (const int i) const
	Get proposed parameter value. More...
double	GetParCurr (const int i) const
	Get current parameter value. More...
const std::vector< double > &	GetParCurrVec () const
	Get vector of current parameter values. More...
double	GetParPreFit (const int i) const
	Get prior parameter value. More...
double	GetUpperBound (const int i) const
	Get upper parameter bound in which it is physically valid. More...
double	GetLowerBound (const int i) const
	Get lower parameter bound in which it is physically valid. More...
double	GetIndivStepScale (const int ParameterIndex) const
	Get individual step scale for selected parameter. More...
double	GetGlobalStepScale () const
	Get global step scale for covariance object. More...
int	GetNParameters () const
	Get number of params which will be different depending if using Eigen decomposition or not. More...
void	PrintPreFitValues () const
	Print prior value for every parameter. More...
void	PrintPreFitCurrPropValues () const
	Print prior, current and proposed value for each parameter. More...
void	PrintIndivStepScale () const
	Print step scale for each parameter. More...
virtual void	ProposeStep ()
	Generate a new proposed state. More...
void	Randomize () _noexcept_
	"Randomize" the parameters in the covariance class for the proposed step. Used the proposal kernel and the current parameter value to set proposed step More...
void	CorrelateSteps () _noexcept_
	Use Cholesky throw matrix for better step proposal. More...
void	UpdateAdaptiveCovariance ()
	Method to update adaptive MCMC [15]. More...
void	AcceptStep () _noexcept_
	Accepted this step. More...
void	SetFixAllParameters ()
	Set all parameters to be fixed at prior values. More...
void	SetFixParameter (const int i)
	Set parameter to be fixed at prior value. More...
void	SetFixParameter (const std::string &name)
	Set parameter to be fixed at prior value. More...
void	SetFreeAllParameters ()
	Set all parameters to be treated as free. More...
void	SetFreeParameter (const int i)
	Set parameter to be treated as free. More...
void	SetFreeParameter (const std::string &name)
	Set parameter to be treated as free. More...
bool	IsParameterFixed (const int i) const
	Is parameter fixed or not. More...
bool	IsParameterFixed (const std::string &name) const
	Is parameter fixed or not. More...
void	ConstructPCA (const double eigen_threshold, int FirstPCAdpar, int LastPCAdpar)
	CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold. More...
bool	IsPCA () const
	is PCA, can use to query e.g. LLH scans More...
YAML::Node	GetConfig () const
	Getter to return a copy of the YAML node. More...
AdaptiveMCMCHandler *	GetAdaptiveHandler () const
	Get pointer for AdaptiveHandler. More...
void	SetTune (const std::string &TuneName)
	KS: Set proposed parameter values vector to be base on tune values, for example set proposed values to be of generated or maybe PostND. More...
PCAHandler *	GetPCAHandler () const
	Get pointer for PCAHandler. More...
void	MatchMaCh3OutputBranches (TTree *PosteriorFile, std::vector< double > &BranchValues, std::vector< std::string > &BranchNames, const std::vector< std::string > &FancyNames={})
	Matches branches in a TTree to parameters in a systematic handler. More...

Protected Member Functions
void	InitFromFile (const std::string &name, const std::string &file)
	Initialisation of the class using matrix from root file. More...
void	ReserveMemory (const int size)
	Initialise vectors with parameters information. More...
void	MakePosDef (TMatrixDSym *cov=nullptr, bool verbose=true)
	Make matrix positive definite by adding small values to diagonal, necessary for inverting matrix. More...
void	MakeClosestPosDef (TMatrixDSym *cov)
	HW: Finds closest possible positive definite matrix in Frobenius Norm ||.||_frob Where ||X||_frob=sqrt[sum_ij(x_ij^2)] (basically just turns an n,n matrix into vector in n^2 space then does Euclidean norm) More...
void	SetThrowMatrixFromFile (const std::string &matrix_file_name, const std::string &matrix_name, const std::string &means_name)
	sets throw matrix from a file More...
void	SanitizeAdaption () const
	Perform sanity check to ensure adaption isn't misbehaving before fit starts. More...
void	FlipParameterValue (const int index, const double FlipPoint)
	KS: Flip parameter around given value, for example mass ordering around 0. More...
void	CircularParBounds (const int i, const double LowBound, const double UpBound)
	HW :: This method is a tad hacky but modular arithmetic gives me a headache. More...
void	EnableSpecialProposal (const YAML::Node &param, const int Index)
	Enable special proposal. More...
void	SpecialStepProposal ()
	Perform Special Step Proposal. More...

Protected Attributes
bool	doSpecialStepProposal
	Check if any of special step proposal were enabled. More...
std::string	inputFile
	The input root file we read in. More...
std::string	matrixName
	Name of cov matrix. More...
TMatrixDSym *	covMatrix
	The covariance matrix. More...
TMatrixDSym *	invCovMatrix
	The inverse covariance matrix. More...
std::vector< std::vector< double > >	InvertCovMatrix
	KS: Same as above but much faster as TMatrixDSym cache miss. More...
std::vector< std::unique_ptr< TRandom3 > >	random_number
	KS: Set Random numbers for each thread so each thread has different seed. More...
double *	randParams
	Random number taken from gaussian around prior error used for corr_throw. More...
double *	corr_throw
	Result of multiplication of Cholesky matrix and randParams. More...
double	_fGlobalStepScale
	Global step scale applied to all params in this class. More...
int	PrintLength
	KS: This is used when printing parameters, sometimes we have super long parameters name, we want to flexibly adjust couts. More...
std::vector< std::string >	_fNames
	ETA _fNames is set automatically in the covariance class to be something like param_i, this is currently to make things compatible with the Diagnostic tools. More...
std::vector< std::string >	_fFancyNames
	Fancy name for example rather than param_0 it is MAQE, useful for human reading. More...
YAML::Node	_fYAMLDoc
	Stores config describing systematics. More...
int	_fNumPar
	Number of systematic parameters. More...
std::vector< double >	_fPreFitValue
	Parameter value dictated by the prior model. Based on it penalty term is calculated. More...
std::vector< double >	_fCurrVal
	Current value of the parameter. More...
std::vector< M3::float_t >	_fPropVal
	Proposed value of the parameter. More...
std::vector< double >	_fError
	Prior error on the parameter. More...
std::vector< double >	_fLowBound
	Lowest physical bound, parameter will not be able to go beyond it. More...
std::vector< double >	_fUpBound
	Upper physical bound, parameter will not be able to go beyond it. More...
std::vector< double >	_fIndivStepScale
	Individual step scale used by MCMC algorithm. More...
std::vector< bool >	_fFlatPrior
	Whether to apply flat prior or not. More...
std::vector< double >	_fIndivStepScaleInitial
	Backup of _fIndivStepScale for parameters which are skipped during adaption. More...
double	_fGlobalStepScaleInitial
	Backup of _fGlobalStepScale for parameters which are skipped during adaption. More...
std::vector< bool >	param_skip_adapt_flags
	Flags telling if parameter should be skipped during adaption. More...
TMatrixDSym *	throwMatrix
	Matrix which we use for step proposal before Cholesky decomposition (not actually used for step proposal) More...
double **	throwMatrixCholDecomp
	Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss. More...
bool	pca
	perform PCA or not More...
bool	use_adaptive
	Are we using AMCMC? More...
std::unique_ptr< PCAHandler >	PCAObj
	Struct containing information about PCA. More...
std::unique_ptr< AdaptiveMCMCHandler >	AdaptiveHandler
	Struct containing information about adaption. More...
std::unique_ptr< ParameterTunes >	Tunes
	Struct containing information about adaption. More...
std::vector< int >	FlipParameterIndex
	Indices of parameters with flip symmetry. More...
std::vector< double >	FlipParameterPoint
	Central points around which parameters are flipped. More...
std::vector< int >	CircularBoundsIndex
	Indices of parameters with circular bounds. More...
std::vector< std::pair< double, double > >	CircularBoundsValues
	Circular bounds for each parameter (lower, upper) More...

Private Member Functions
void	ToggleFixParameter (const int i)
	Toggle fixing parameter at prior values. More...
void	ToggleFixParameter (const std::string &name)
	Toggle fixing parameter at prior values. More...

Detailed Description

Base class for handling systematic uncertainty parameters.

Provides core functionality for managing systematic parameters, including likelihood evaluation and covariance handling.

How parameters are loaded (e.g. from configuration files or ROOT inputs) is left to derived classes.

In the context of MCMC, the class may also be responsible for proposing parameter steps. Proposal strategies may include PCA-based sampling or adaptive proposal matrices.

Author

Dan Barrow

Ed Atkin

Kamil Skwarczynski

Definition at line 25 of file ParameterHandlerBase.h.

Constructor & Destructor Documentation

ParameterHandlerBase() [1/2]

ParameterHandlerBase::ParameterHandlerBase	(	std::string	name,
		std::string	file,
		double	threshold = -1,
		int	FirstPCAdpar = -999,
		int	LastPCAdpar = -999
	)

"Usual" constructors from root file

Parameters

name	Matrix name
file	Path to matrix root file

Definition at line 7 of file ParameterHandlerBase.cpp.

                     : inputFile(file), pca(true) {
// ********************************************
  MACH3LOG_DEBUG("Constructing instance of ParameterHandler");
  doSpecialStepProposal = false;
  // Not using adaptive by default
  use_adaptive = false;
  if (threshold < 0 || threshold >= 1) {
    MACH3LOG_INFO("NOTE: {} {}", name, file);
    MACH3LOG_INFO("Principal component analysis but given the threshold for the principal components to be less than 0, or greater than (or equal to) 1. This will not work");
    MACH3LOG_INFO("Please specify a number between 0 and 1");
    MACH3LOG_INFO("You specified: ");
    MACH3LOG_INFO("Am instead calling the usual non-PCA constructor...");
    pca = false;
  }
 
  InitFromFile(name, file);
 
  // Call the innocent helper function
  if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
}

ParameterHandlerBase() [2/2]

ParameterHandlerBase::ParameterHandlerBase ( )

default

~ParameterHandlerBase()

ParameterHandlerBase::~ParameterHandlerBase ( )

virtual

Destructor.

Definition at line 31 of file ParameterHandlerBase.cpp.

                                           {
// ********************************************
  delete[] randParams;
  delete[] corr_throw;
 
  if (covMatrix != nullptr) delete covMatrix;
  if (invCovMatrix != nullptr) delete invCovMatrix;
  if (throwMatrix != nullptr) delete throwMatrix;
  for(int i = 0; i < _fNumPar; i++) {
    delete[] throwMatrixCholDecomp[i];
  }
  delete[] throwMatrixCholDecomp;
}

Member Function Documentation

AcceptStep()

void ParameterHandlerBase::AcceptStep

(

)

Accepted this step.

Definition at line 503 of file ParameterHandlerBase.cpp.

                                                 {
// ********************************************
  if (!pca) {
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < _fNumPar; ++i) {
      // Update state so that current state is proposed state
      _fCurrVal[i] = _fPropVal[i];
    }
  } else {
    PCAObj->AcceptStep();
  }
 
  if (AdaptiveHandler) {
    AdaptiveHandler->IncrementAcceptedSteps();
  }
}

CalcLikelihood()

double ParameterHandlerBase::CalcLikelihood

(

)

const

Calc penalty term based on inverted covariance matrix.

The log-likelihood is computed as:

\[ \log \mathcal{L} = \frac{1}{2} \sum_{i}^{\textrm{pars}} \sum_{j}^{\textrm{pars}} \Delta \vec{p}_i \left( V^{-1} \right)_{i,j} \Delta \vec{p}_j \]

where:

• \(\Delta \vec{p}_i = \theta_i - \theta_{i,0}\) is the difference between the current and pre-fit parameter values,
• \(V^{-1}\) is the inverted covariance matrix.

Note

• If _fFlatPrior[i] is true, the parameter is excluded from the calculation.

Definition at line 573 of file ParameterHandlerBase.cpp.

                                                             {
// ********************************************
  double logL = 0.0;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:logL)
  #endif
  for(int i = 0; i < _fNumPar; ++i) {
    if(_fFlatPrior[i]){
      //HW: Flat prior, no need to calculate anything
      continue;
    }
    // KS: Precalculate Diff once per "i" without doing this for every "j"
    const double Diff = _fPropVal[i] - _fPreFitValue[i];
    #ifdef MULTITHREAD
    #pragma omp simd
    #endif
    for (int j = 0; j <= i; ++j) {
      if (!_fFlatPrior[j]) {
        //KS: Since matrix is symmetric we can calculate non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
        double scale = (i != j) ? 1. : 0.5;
        logL += scale * Diff * (_fPropVal[j] - _fPreFitValue[j])*InvertCovMatrix[i][j];
      }
    }
  }
  return logL;
}

CheckBounds()

int ParameterHandlerBase::CheckBounds

(

)

const

Check if parameters were proposed outside physical boundary.

Definition at line 601 of file ParameterHandlerBase.cpp.

                                                       {
// ********************************************
  int NOutside = 0;
  for (int i = 0; i < _fNumPar; ++i) {
    // KS: Count how many parameters are outside bounds using branchless logic
    // faster by at least factor two
    // Do not multithread even with 5k params no gains
    NOutside += (_fPropVal[i] > _fUpBound[i]) | (_fPropVal[i] < _fLowBound[i]);
  }
  return NOutside;
}

CircularParBounds()

void ParameterHandlerBase::CircularParBounds ( const int  i, const double  LowBound, const double  UpBound  )

protected

HW :: This method is a tad hacky but modular arithmetic gives me a headache.

Author

Henry Wallace

Definition at line 526 of file ParameterHandlerBase.cpp.

                                                                                                         {
// *************************************
  if(_fPropVal[index] > UpBound) {
    _fPropVal[index] = static_cast<M3::float_t>(LowBound + std::fmod(_fPropVal[index] - UpBound, UpBound - LowBound));
  } else if (_fPropVal[index] < LowBound) {
    _fPropVal[index] = static_cast<M3::float_t>(UpBound - std::fmod(LowBound - _fPropVal[index], UpBound - LowBound));
  }
}

ConstructPCA()

void ParameterHandlerBase::ConstructPCA	(	const double	eigen_threshold,
		int	FirstPCAdpar,
		int	LastPCAdpar
	)

CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold.

Parameters

eigen_threshold	PCA threshold from 0 to 1. Default is -1 and means no PCA
FirstPCAdpar	First PCA parameter that will be decomposed.
LastPCAdpar	First PCA parameter that will be decomposed.

See also

For more details, visit the PCAHandler

Definition at line 46 of file ParameterHandlerBase.cpp.

                                                                                                       {
// ********************************************
  if(AdaptiveHandler) {
    MACH3LOG_ERROR("Adaption has been enabled and now trying to enable PCA. Right now both configuration don't work with each other");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  PCAObj = std::make_unique<PCAHandler>();
  //Check whether first and last pcadpar are set and if not just PCA everything
  if(FirstPCAdpar == -999 || LastPCAdpar == -999) {
    if(FirstPCAdpar == -999 && LastPCAdpar == -999) {
      FirstPCAdpar = 0;
      LastPCAdpar = covMatrix->GetNrows()-1;
    }
    else{
      MACH3LOG_ERROR("You must either leave FirstPCAdpar and LastPCAdpar at -999 or set them both to something");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  }
 
  PCAObj->ConstructPCA(covMatrix, FirstPCAdpar, LastPCAdpar, eigen_threshold, _fNumPar);
  PCAObj->SetupPointers(&_fCurrVal, &_fPropVal);
  // Make a note that we have now done PCA
  pca = true;
}

CorrelateSteps()

void ParameterHandlerBase::CorrelateSteps

(

)

Use Cholesky throw matrix for better step proposal.

Definition at line 478 of file ParameterHandlerBase.cpp.

                                                     {
// ************************************************
  //KS: Using custom function compared to ROOT one with 8 threads we have almost factor 2 performance increase, by replacing TMatrix with just double we increase it even more
  M3::MatrixVectorMulti(corr_throw, throwMatrixCholDecomp, randParams, _fNumPar);
 
  // If not doing PCA
  if (!pca) {
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < _fNumPar; ++i) {
      if (!IsParameterFixed(i) > 0.) {
        #pragma GCC diagnostic push
        #pragma GCC diagnostic ignored "-Wuseless-cast"
        _fPropVal[i] = static_cast<M3::float_t>(_fCurrVal[i] + corr_throw[i]*_fGlobalStepScale*_fIndivStepScale[i]);
        #pragma GCC diagnostic pop
      }
    }
    // If doing PCA throw uncorrelated in PCA basis (orthogonal basis by definition)
  } else {
    PCAObj->CorrelateSteps(_fIndivStepScale, _fGlobalStepScale, randParams, corr_throw);
  }
}

EnableSpecialProposal()

void ParameterHandlerBase::EnableSpecialProposal ( const YAML::Node &  param, const int  Index  )

protected

Enable special proposal.

Definition at line 116 of file ParameterHandlerBase.cpp.

                                                                                      {
// ********************************************
  doSpecialStepProposal = true;
 
  bool CircEnabled = false;
  bool FlipEnabled = false;
 
  if (param["CircularBounds"]) {
    CircEnabled = true;
  }
 
  if (param["FlipParameter"]) {
    FlipEnabled = true;
  }
 
  if (!CircEnabled && !FlipEnabled) {
    MACH3LOG_ERROR("None of Special Proposal were enabled even though param {}, has SpecialProposal entry in Yaml", GetParFancyName(Index));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (CircEnabled) {
    CircularBoundsIndex.push_back(Index);
    CircularBoundsValues.push_back(Get<std::pair<double, double>>(param["CircularBounds"], __FILE__, __LINE__));
    MACH3LOG_INFO("Enabling CircularBounds for parameter {} with range [{}, {}]",
                  GetParFancyName(Index),
                  CircularBoundsValues.back().first,
                  CircularBoundsValues.back().second);
    // KS: Make sure circular bounds are within physical bounds. If we are outside of physics bound MCMC will never explore such phase space region
    if (CircularBoundsValues.back().first < _fLowBound.at(Index) || CircularBoundsValues.back().second > _fUpBound.at(Index)) {
      MACH3LOG_ERROR("Circular bounds [{}, {}] for parameter {} exceed physical bounds [{}, {}]",
                     CircularBoundsValues.back().first, CircularBoundsValues.back().second,
                     GetParFancyName(Index),
                     _fLowBound.at(Index), _fUpBound.at(Index));
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    // KS: Make sure CircularPrior is applied only to param with flat prior. Sadly doesn't work with Gaussian
    if(GetFlatPrior(Index) == false) {
      MACH3LOG_ERROR("Enabled CircularPrior for parameter {}, which has gaussian prior", GetParFancyName(Index));
      MACH3LOG_ERROR("This is not supported, CircularPrior only works with flat prior");
      MACH3LOG_ERROR("Change FlatPrior in Parameter config to true");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  if (FlipEnabled) {
    FlipParameterIndex.push_back(Index);
    FlipParameterPoint.push_back(Get<double>(param["FlipParameter"], __FILE__, __LINE__));
    MACH3LOG_INFO("Enabling Flipping for parameter {} with value {}",
                  GetParFancyName(Index),
                  FlipParameterPoint.back());
  }
 
  if (CircEnabled && FlipEnabled) {
    if (FlipParameterPoint.back() < CircularBoundsValues.back().first || FlipParameterPoint.back() > CircularBoundsValues.back().second) {
      MACH3LOG_ERROR("FlipParameter value {} for parameter {} is outside the CircularBounds [{}, {}]",
                     FlipParameterPoint.back(), GetParFancyName(Index), CircularBoundsValues.back().first, CircularBoundsValues.back().second);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
 
    const double low = CircularBoundsValues.back().first;
    const double high = CircularBoundsValues.back().second;
 
    // Sanity check: ensure flipping any x in [low, high] keeps the result in [low, high]
    const double flipped_low = 2 * FlipParameterPoint.back() - low;
    const double flipped_high = 2 * FlipParameterPoint.back() - high;
    const double min_flip = std::min(flipped_low, flipped_high);
    const double max_flip = std::max(flipped_low, flipped_high);
 
    if (min_flip < low || max_flip > high) {
      MACH3LOG_ERROR("Flipping about point {} for parameter {} would leave circular bounds [{}, {}]",
                     FlipParameterPoint.back(), GetParFancyName(Index), low, high);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
}

FlipParameterValue()

void ParameterHandlerBase::FlipParameterValue ( const int  index, const double  FlipPoint  )

protected

KS: Flip parameter around given value, for example mass ordering around 0.

Parameters

index	parameter index you want to flip
FlipPoint	Value around which flipping is done

Definition at line 536 of file ParameterHandlerBase.cpp.

                                                                                     {
// *************************************
  if(random_number[0]->Uniform() < 0.5) {
    _fPropVal[index] = static_cast<M3::float_t>(2 * FlipPoint - _fPropVal[index]);
  }
}

GetAdaptiveHandler()

AdaptiveMCMCHandler* ParameterHandlerBase::GetAdaptiveHandler ( ) const

inline

Get pointer for AdaptiveHandler.

Definition at line 310 of file ParameterHandlerBase.h.

                                                   {
    if (!use_adaptive) {
      MACH3LOG_ERROR("Am not running in Adaptive mode");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    return AdaptiveHandler.get();
  }

GetConfig()

YAML::Node ParameterHandlerBase::GetConfig ( ) const

inline

Getter to return a copy of the YAML node.

Definition at line 307 of file ParameterHandlerBase.h.

{ return _fYAMLDoc; }

GetCorrelationMatrix()

TH2D * ParameterHandlerBase::GetCorrelationMatrix

(

)

const

KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting.

This function converts the covariance matrix to a correlation matrix and returns a TH2D object, which can be used for advanced plotting purposes.

Returns

A pointer to a TH2D object representing the correlation matrix

Definition at line 1168 of file ParameterHandlerBase.cpp.

                                                       {
// ********************************************
  TH2D* hMatrix = new TH2D(GetName().c_str(), GetName().c_str(), _fNumPar, 0.0, _fNumPar, _fNumPar, 0.0, _fNumPar);
  hMatrix->SetDirectory(nullptr);
  for(int i = 0; i < _fNumPar; i++)
  {
    hMatrix->SetBinContent(i+1, i+1, 1.);
    hMatrix->GetXaxis()->SetBinLabel(i+1, GetParFancyName(i).c_str());
    hMatrix->GetYaxis()->SetBinLabel(i+1, GetParFancyName(i).c_str());
  }
 
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for(int i = 0; i < _fNumPar; i++)
  {
    for(int j = 0; j <= i; j++)
    {
      const double Corr = (*covMatrix)(i,j) / ( GetDiagonalError(i) * GetDiagonalError(j));
      hMatrix->SetBinContent(i+1, j+1, Corr);
      hMatrix->SetBinContent(j+1, i+1, Corr);
    }
  }
  return hMatrix;
}

GetCorrThrows()

double ParameterHandlerBase::GetCorrThrows ( const int  i ) const

inline

Return correlated throws.

Parameters

i	Parameter index

Definition at line 134 of file ParameterHandlerBase.h.

{ return corr_throw[i]; }

GetCovMatrix()

TMatrixDSym* ParameterHandlerBase::GetCovMatrix ( ) const

inline

Return covariance matrix.

Definition at line 126 of file ParameterHandlerBase.h.

{ return covMatrix; }

GetDiagonalError()

double ParameterHandlerBase::GetDiagonalError ( const int  i ) const

inline

Get diagonal error for ith parameter.

Parameters

i	Parameter index

Definition at line 155 of file ParameterHandlerBase.h.

{ return std::sqrt((*covMatrix)(i,i)); }

GetDoAdaption()

bool ParameterHandlerBase::GetDoAdaption ( ) const

inline

Do we adapt or not.

Definition at line 174 of file ParameterHandlerBase.h.

{return use_adaptive;}

GetError()

double ParameterHandlerBase::GetError ( const int  i ) const

inline

Get the error for the ith parameter.

Parameters

i	Parameter index

Definition at line 158 of file ParameterHandlerBase.h.

{return _fError[i];}

GetFlatPrior()

bool ParameterHandlerBase::GetFlatPrior ( const int  i ) const

inline

Get if param has flat prior or not.

Parameters

i	Parameter index

Definition at line 138 of file ParameterHandlerBase.h.

{ return _fFlatPrior[i]; }

GetGlobalStepScale()

double ParameterHandlerBase::GetGlobalStepScale ( ) const

inline

Get global step scale for covariance object.

Definition at line 240 of file ParameterHandlerBase.h.

{return _fGlobalStepScale; }

GetIndivStepScale()

double ParameterHandlerBase::GetIndivStepScale ( const int  ParameterIndex ) const

inline

Get individual step scale for selected parameter.

Parameters

ParameterIndex

Parameter index

Definition at line 238 of file ParameterHandlerBase.h.

{return _fIndivStepScale.at(ParameterIndex); }

GetInvCovMatrix() [1/2]

TMatrixDSym* ParameterHandlerBase::GetInvCovMatrix ( ) const

inline

Return inverted covariance matrix.

Definition at line 128 of file ParameterHandlerBase.h.

{ return invCovMatrix; }

GetInvCovMatrix() [2/2]

double ParameterHandlerBase::GetInvCovMatrix ( const int  i, const int  j  ) const

inline

Return inverted covariance matrix.

Definition at line 130 of file ParameterHandlerBase.h.

{ return InvertCovMatrix[i][j]; }

GetLikelihood()

double ParameterHandlerBase::GetLikelihood ( )

virtual

Return CalcLikelihood if some params were thrown out of boundary return LARGE_LOGL

Reimplemented in PyParameterHandlerBase.

Definition at line 614 of file ParameterHandlerBase.cpp.

                                           {
// ********************************************
  // Default behaviour is to reject negative values + do std llh calculation
  const int NOutside = CheckBounds();
 
  if(NOutside > 0) return NOutside*M3::_LARGE_LOGL_;
 
  return CalcLikelihood();
}

GetLowerBound()

double ParameterHandlerBase::GetLowerBound ( const int  i ) const

inline

Get lower parameter bound in which it is physically valid.

Parameters

i	Parameter index

Definition at line 235 of file ParameterHandlerBase.h.

{ return _fLowBound[i]; }

GetName()

std::string ParameterHandlerBase::GetName ( ) const

inline

Get name of covariance.

Definition at line 141 of file ParameterHandlerBase.h.

{ return matrixName; }

GetNParameters()

int ParameterHandlerBase::GetNParameters ( ) const

inline

Get number of params which will be different depending if using Eigen decomposition or not.

Definition at line 243 of file ParameterHandlerBase.h.

                             {
    if (pca) return PCAObj->GetNumberPCAedParameters();
    else return _fNumPar;
  }

GetNumParams()

int ParameterHandlerBase::GetNumParams ( ) const

inline

Get total number of parameters.

Definition at line 213 of file ParameterHandlerBase.h.

{return _fNumPar;}

GetParCurr()

double ParameterHandlerBase::GetParCurr ( const int  i ) const

inline

Get current parameter value.

Parameters

i	Parameter index

Definition at line 223 of file ParameterHandlerBase.h.

{ return _fCurrVal[i]; }

GetParCurrVec()

const std::vector<double>& ParameterHandlerBase::GetParCurrVec ( ) const

inline

Get vector of current parameter values.

Definition at line 225 of file ParameterHandlerBase.h.

{ return _fCurrVal; }

GetParFancyName()

std::string ParameterHandlerBase::GetParFancyName ( const int  i ) const

inline

Get fancy name of the Parameter.

Parameters

i	Parameter index

Definition at line 151 of file ParameterHandlerBase.h.

{return _fFancyNames[i];}

GetParIndex()

int ParameterHandlerBase::GetParIndex

(

const std::string &

name

)

const

Get index based on name.

Definition at line 703 of file ParameterHandlerBase.cpp.

                                                                 {
// ********************************************
  int Index = M3::_BAD_INT_;
  for (int i = 0; i <_fNumPar; ++i) {
    if(name == _fFancyNames[i]) {
      Index = i;
      break;
    }
  }
  return Index;
}

GetParName()

std::string ParameterHandlerBase::GetParName ( const int  i ) const

inline

Get name of parameter.

Parameters

i	Parameter index

Definition at line 144 of file ParameterHandlerBase.h.

{return _fNames[i];}

GetParPreFit()

double ParameterHandlerBase::GetParPreFit ( const int  i ) const

inline

Get prior parameter value.

Parameters

i	Parameter index

Definition at line 229 of file ParameterHandlerBase.h.

{ return _fPreFitValue[i]; }

GetParProp()

M3::float_t ParameterHandlerBase::GetParProp ( const int  i ) const

inline

Get proposed parameter value.

Parameters

i	Parameter index

Definition at line 220 of file ParameterHandlerBase.h.

{ return _fPropVal[i]; }

GetParPropVec()

const std::vector<M3::float_t>& ParameterHandlerBase::GetParPropVec ( )

inline

Get a reference to the proposed parameter values Can be useful if you want to track these without having to copy values using getProposed()

Definition at line 210 of file ParameterHandlerBase.h.

{return _fPropVal;}

GetPCAHandler()

PCAHandler* ParameterHandlerBase::GetPCAHandler ( ) const

inline

Get pointer for PCAHandler.

Definition at line 322 of file ParameterHandlerBase.h.

                                    {
    if (!pca) {
      MACH3LOG_ERROR("Am not running in PCA mode");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    return PCAObj.get();
  }

GetPreFitValues()

std::vector<double> ParameterHandlerBase::GetPreFitValues ( ) const

inline

Get the pre-fit values of the parameters.

Definition at line 215 of file ParameterHandlerBase.h.

{return _fPreFitValue;}

GetProposed()

std::vector< double > ParameterHandlerBase::GetProposed

(

)

const

Get vector of all proposed parameter values.

Definition at line 275 of file ParameterHandlerBase.cpp.

                                                          {
// ********************************************
  std::vector<double> props(_fNumPar);
  for (int i = 0; i < _fNumPar; ++i) props[i] = _fPropVal[i];
  return props;
}

GetRandomThrow()

double ParameterHandlerBase::GetRandomThrow ( const int  i ) const

inline

Get random value useful for debugging/CI.

Parameters

i	Parameter index

Definition at line 77 of file ParameterHandlerBase.h.

{ return randParams[i];}

GetThrowMatrix() [1/2]

TMatrixDSym* ParameterHandlerBase::GetThrowMatrix ( ) const

inline

Get matrix used for step proposal.

Definition at line 187 of file ParameterHandlerBase.h.

{return throwMatrix;}

GetThrowMatrix() [2/2]

double ParameterHandlerBase::GetThrowMatrix ( const int  i, const int  j  ) const

inline

Get matrix used for step proposal.

Definition at line 189 of file ParameterHandlerBase.h.

{ return throwMatrixCholDecomp[i][j];}

GetUpperBound()

double ParameterHandlerBase::GetUpperBound ( const int  i ) const

inline

Get upper parameter bound in which it is physically valid.

Parameters

i	Parameter index

Definition at line 232 of file ParameterHandlerBase.h.

{ return _fUpBound[i];}

InitFromFile()

void ParameterHandlerBase::InitFromFile ( const std::string &  name, const std::string &  file  )

protected

Initialisation of the class using matrix from root file.

Definition at line 73 of file ParameterHandlerBase.cpp.

                                                                                    {
// ********************************************
  // Set the covariance matrix from input ROOT file (e.g. flux, ND280, NIWG)
  TFile *infile = new TFile(file.c_str(), "READ");
  if (infile->IsZombie()) {
    MACH3LOG_ERROR("Could not open input covariance ROOT file {} !!!", file);
    MACH3LOG_ERROR("Was about to retrieve matrix with name {}", name);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  TMatrixDSym *CovMat = static_cast<TMatrixDSym*>(infile->Get(name.c_str()));
 
  if (!CovMat) {
    MACH3LOG_ERROR("Could not find covariance matrix name {} in file {}", name, file);
    MACH3LOG_ERROR("Are you really sure {} exists in the file?", name);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  PrintLength = 35;
 
  const int nThreads = M3::GetNThreads();
  //KS: set Random numbers for each thread so each thread has different seed
  //or for one thread if without MULTITHREAD
  random_number.reserve(nThreads);
  for (int iThread = 0; iThread < nThreads; iThread++) {
    random_number.emplace_back(std::make_unique<TRandom3>(0));
  }
  // Set the covariance matrix
  _fNumPar = CovMat->GetNrows();
 
  ReserveMemory(_fNumPar);
  SetName(name);
  MakePosDef(CovMat);
  SetCovMatrix(CovMat);
 
  infile->Close();
 
  MACH3LOG_INFO("Created covariance matrix named: {}", GetName());
  MACH3LOG_INFO("from file: {}", file);
  delete infile;
}

InitialiseAdaption()

void ParameterHandlerBase::InitialiseAdaption

(

const YAML::Node &

adapt_manager

)

Initialise adaptive MCMC.

Parameters

adapt_manager

Node having from which we load all adaptation options

Definition at line 999 of file ParameterHandlerBase.cpp.

                                                                           {
// ********************************************
  if(PCAObj){
    MACH3LOG_ERROR("PCA has been enabled and now trying to enable Adaption. Right now both configuration don't work with each other");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  if(AdaptiveHandler){
    MACH3LOG_ERROR("Adaptive Handler has already been initialise can't do it again so skipping.");
    return;
  }
  AdaptiveHandler = std::make_unique<AdaptiveMCMCHandler>();
 
  // HH: Backing up _fIndivStepScale and _fGlobalStepScale before adaption
  _fIndivStepScaleInitial = _fIndivStepScale;
  _fGlobalStepScaleInitial = _fGlobalStepScale;
 
  // HH: adding these here because they will be used to set the individual step scales for non-adapting parameters
  auto params_to_skip = GetFromManager<std::vector<std::string>>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ParametersToSkip"], {}, __FILE__ , __LINE__);
  // Build a list of skip flags
  param_skip_adapt_flags.resize(_fNumPar, false);
  for (int i = 0; i <_fNumPar; ++i) {
    param_skip_adapt_flags[i] = M3::CaseInsensitiveMatchAny(_fFancyNames[i], params_to_skip);
  }
 
  // Now we read the general settings [these SHOULD be common across all matrices!]
  bool success = AdaptiveHandler->InitFromConfig(adapt_manager, matrixName,
    &_fFancyNames, &_fCurrVal, &_fError,
    &param_skip_adapt_flags, throwMatrix, _fGlobalStepScaleInitial
  );
  if (success) {
    // Ensure there is no misconfiguration in adaption config
    SanitizeAdaption();
    AdaptiveHandler->Print();
  } else {
    MACH3LOG_INFO("Not using adaptive MCMC for {}. Checking external matrix options...", matrixName);
  }
 
  // HH: Adjusting the external matrix reading logic such that you can not do adaptive
  // and still read an external matrix
  // Logic:
  // if read external matrix:
  //    set throw matrix regardless of adaptive or not
  // else:
  //    if adaptive:
  //       create new adaptive matrix from scratch
  //    else:
  //       do nothing
  if(GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["UseExternalMatrix"], false, __FILE__ , __LINE__)) {
    // Finally, we accept that we want to read the matrix from a file!
    auto external_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixFileName"], "", __FILE__ , __LINE__);
    auto external_matrix_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixName"], "", __FILE__ , __LINE__);
    auto external_mean_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMeansName"], "", __FILE__ , __LINE__);
 
    AdaptiveHandler->SetThrowMatrixFromFile(external_file_name, external_matrix_name, external_mean_name, use_adaptive);
    SetThrowMatrix(AdaptiveHandler->GetAdaptiveCovariance());
 
    ResetIndivStepScale();
    // HH: Set individual step scales for non-adapting parameters to the default individual step scales
    // global step scale should be 1 so no need to adjust for that
    SetIndivStepScaleForSkippedAdaptParams();
 
    MACH3LOG_INFO("Successfully Set External Throw Matrix Stored in {}", external_file_name);
  } else {
    MACH3LOG_INFO("Not using external matrix for {}", matrixName);
    if (!success) return; // Not adaptive either so nothing to do
    MACH3LOG_INFO("Initialising adaption from scratch");
    // If we don't have a covariance matrix to start from for adaptive tune we need to make one!
    use_adaptive = true;
    AdaptiveHandler->CheckMatrixValidityForAdaption(GetCovMatrix());
    AdaptiveHandler->CreateNewAdaptiveCovariance();
    return;
  }
}

IsParameterFixed() [1/2]

bool ParameterHandlerBase::IsParameterFixed ( const int  i ) const

inline

Is parameter fixed or not.

Parameters

i	Parameter index

Definition at line 288 of file ParameterHandlerBase.h.

                                           {
    if (_fError[i] < 0) { return true; }
    else                { return false; }
  }

IsParameterFixed() [2/2]

bool ParameterHandlerBase::IsParameterFixed

(

const std::string &

name

)

const

Is parameter fixed or not.

Parameters

name	Name of parameter you want to check if is fixed

Definition at line 793 of file ParameterHandlerBase.cpp.

                                                                       {
// ********************************************
  const int Index = GetParIndex(name);
  if(Index != M3::_BAD_INT_) {
    return IsParameterFixed(Index);
  }
 
  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore don't know if it fixed", name);
  return false;
}

IsPCA()

bool ParameterHandlerBase::IsPCA ( ) const

inline

is PCA, can use to query e.g. LLH scans

Definition at line 304 of file ParameterHandlerBase.h.

{ return pca; }

MakeClosestPosDef()

void ParameterHandlerBase::MakeClosestPosDef ( TMatrixDSym *  cov )

protected

HW: Finds closest possible positive definite matrix in Frobenius Norm ||.||_frob Where ||X||_frob=sqrt[sum_ij(x_ij^2)] (basically just turns an n,n matrix into vector in n^2 space then does Euclidean norm)

Definition at line 1121 of file ParameterHandlerBase.cpp.

                                                             {
// ********************************************
  // Want to get cov' = (cov_sym+cov_polar)/2
  // cov_sym=(cov+cov^T)/2
  // cov_polar-> SVD cov to cov=USV^T then cov_polar=VSV^T
 
  //Get frob norm of cov
  //  Double_t cov_norm=cov->E2Norm();
 
  TMatrixDSym* cov_trans = cov;
  cov_trans->T();
  TMatrixDSym cov_sym = 0.5*(*cov+*cov_trans); //If cov is symmetric does nothing, otherwise just ensures symmetry
 
  //Do SVD to get polar form
  TDecompSVD cov_sym_svd=TDecompSVD(cov_sym);
  if(!cov_sym_svd.Decompose()){
    MACH3LOG_WARN("Cannot do SVD on input matrix, trying MakePosDef() first!");
    MakePosDef(&cov_sym);
  }
 
  TMatrixD cov_sym_v = cov_sym_svd.GetV();
  TMatrixD cov_sym_vt = cov_sym_v;
  cov_sym_vt.T();
  //SVD returns as vector (grrr) so need to get into matrix form for multiplying!
  TVectorD cov_sym_sigvect = cov_sym_svd.GetSig();
 
  const Int_t nCols = cov_sym_v.GetNcols(); //square so only need rows hence lack of cols
  TMatrixDSym cov_sym_sig(nCols);
  TMatrixDDiag cov_sym_sig_diag(cov_sym_sig);
  cov_sym_sig_diag=cov_sym_sigvect;
 
  //Can finally get H=VSV
  TMatrixDSym cov_sym_polar = cov_sym_sig.SimilarityT(cov_sym_vt);//V*S*V^T (this took forver to find!)
 
  //Now we can construct closest approximater Ahat=0.5*(B+H)
  TMatrixDSym cov_closest_approx  = 0.5*(cov_sym+cov_sym_polar);//Not fully sure why this is even needed since symmetric B -> U=V
  //Get norm of transformed
  //  Double_t approx_norm=cov_closest_approx.E2Norm();
  //MACH3LOG_INFO("Initial Norm: {:.6f} | Norm after transformation: {:.6f} | Ratio: {:.6f}", cov_norm, approx_norm, cov_norm / approx_norm);
 
  *cov = cov_closest_approx;
  //Now can just add a makeposdef!
  MakePosDef(cov);
}

MakePosDef()

void ParameterHandlerBase::MakePosDef ( TMatrixDSym *  cov = nullptr, bool  verbose = true  )

protected

Make matrix positive definite by adding small values to diagonal, necessary for inverting matrix.

Parameters

cov	Matrix which we evaluate Positive Definitiveness
verbose	Do we want it to print warning if it was not positive definite and we had to fix it

Definition at line 853 of file ParameterHandlerBase.cpp.

                                                                    {
// ********************************************
  if(cov == nullptr){
    cov = &*covMatrix;
    MACH3LOG_WARN("Passed nullptr to cov matrix in {}", matrixName);
  }
 
  int n_attempts = M3::MakeMatrixPosDef(cov);
 
  if(n_attempts > 0 && verbose) {
    MACH3LOG_WARN("Covariance matrix {} was not positive-definite, made it positive-definite after {} attempts", matrixName, n_attempts);
  }
 
}

MatchMaCh3OutputBranches()

void ParameterHandlerBase::MatchMaCh3OutputBranches	(	TTree *	PosteriorFile,
		std::vector< double > &	BranchValues,
		std::vector< std::string > &	BranchNames,
		const std::vector< std::string > &	FancyNames = {}
	)

Matches branches in a TTree to parameters in a systematic handler.

Parameters

PosteriorFile	Pointer to the ROOT TTree from MaCh3 fit.
BranchValues	Vector to store the values of the branches (resized inside).
BranchNames	Vector to store the names of the branches (resized inside).
FancyNames	Optional vector of "fancy names" to match. If empty, all parameters are matched.

• If FancyNames is provided, it matches only the parameters whose "fancy names" are in FancyNames. This is useful for studies where one performs ND fits and passes them to FD fits, which may have additional parameters (e.g., oscillations).
• If FancyNames is empty, it matches all parameters in the systematic handler.

Exceptions

MaCh3Exception

if any parameter branch is uninitialized.

Definition at line 1234 of file ParameterHandlerBase.cpp.

                                                                      {
// *************************************
  BranchValues.resize(GetNumParams());
  BranchNames.resize(GetNumParams());
 
  // if fancy names are passed match ONLY them
  // this allow to perform studies where one perform ND fits and pass it to FD fits
  // which usually have more params like osc...
  if(FancyNames.size() != 0){
    for (int i = 0; i < GetNumParams(); ++i) {
      BranchNames[i] = GetParName(i);
      // by default set current step
      BranchValues[i] = _fPropVal[i];
      bool matched = false;
      for (size_t iPar = 0; iPar < FancyNames.size(); ++iPar) {
        if(GetParFancyName(i) == FancyNames[iPar]) {
          MACH3LOG_DEBUG("Matched name {} in config", FancyNames[iPar]);
          PosteriorFile->SetBranchStatus(BranchNames[i].c_str(), true);
          PosteriorFile->SetBranchAddress(BranchNames[i].c_str(), &BranchValues[i]);
          matched = true;
          break;
        }
      }
      if(!matched) {
        MACH3LOG_WARN("Didn't match param {} is this what you want?", GetParFancyName(i));
      }
    }
  } else {
    // simply loop over params and match them
    for (int i = 0; i < GetNumParams(); ++i) {
      BranchNames[i] = GetParName(i);
      if (!PosteriorFile->GetBranch(BranchNames[i].c_str())) {
        MACH3LOG_ERROR("Branch '{}' does not exist in the TTree!", BranchNames[i]);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      PosteriorFile->SetBranchStatus(BranchNames[i].c_str(), true);
      PosteriorFile->SetBranchAddress(BranchNames[i].c_str(), &BranchValues[i]);
    }
  }
}

PrintIndivStepScale()

void ParameterHandlerBase::PrintIndivStepScale

(

)

const

Print step scale for each parameter.

Definition at line 841 of file ParameterHandlerBase.cpp.

                                                     {
// ********************************************
  MACH3LOG_INFO("============================================================");
  MACH3LOG_INFO("{:<{}} | {:<11}", "Parameter:", PrintLength, "Step scale:");
  for (int iParam = 0; iParam < _fNumPar; iParam++) {
    MACH3LOG_INFO("{:<{}} | {:<11}", _fFancyNames[iParam].c_str(), PrintLength, _fIndivStepScale[iParam]);
  }
  MACH3LOG_INFO("============================================================");
}

PrintPreFitCurrPropValues()

void ParameterHandlerBase::PrintPreFitCurrPropValues

(

)

const

Print prior, current and proposed value for each parameter.

Definition at line 555 of file ParameterHandlerBase.cpp.

                                                           {
// ********************************************
  MACH3LOG_INFO("Printing parameters for {}", GetName());
  // Dump out the PCA parameters too
  if (pca) {
    PCAObj->Print();
  }
  MACH3LOG_INFO("{:<30} {:<10} {:<10} {:<10}", "Name", "Prior", "Current", "Proposed");
  for (int i = 0; i < _fNumPar; ++i) {
    MACH3LOG_INFO("{:<30} {:<10.2f} {:<10.2f} {:<10.2f}", GetParFancyName(i), _fPreFitValue[i], _fCurrVal[i], _fPropVal[i]);
  }
}

PrintPreFitValues()

void ParameterHandlerBase::PrintPreFitValues

(

)

const

Print prior value for every parameter.

Definition at line 545 of file ParameterHandlerBase.cpp.

                                                   {
// ********************************************
  MACH3LOG_INFO("Prior values for {} ParameterHandler:", GetName());
  for (int i = 0; i < _fNumPar; i++) {
    MACH3LOG_INFO("    {}   {} ", GetParFancyName(i), GetParPreFit(i));
  }
}

ProposeStep()

void ParameterHandlerBase::ProposeStep ( )

virtual

Generate a new proposed state.

Reimplemented in PyParameterHandlerBase.

Definition at line 404 of file ParameterHandlerBase.cpp.

                                       {
// ************************************************
  // Make the random numbers for the step proposal
  Randomize();
  CorrelateSteps();
 
  // KS: According to Dr Wallace we update using previous not proposed step
  // this way we do special proposal after adaptive after.
  // This way we can shortcut and skip rest of proposal
  if(!doSpecialStepProposal) return;
 
  SpecialStepProposal();
}

RandomConfiguration()

void ParameterHandlerBase::RandomConfiguration

(

)

Randomly throw the parameters in their 1 sigma range.

Definition at line 342 of file ParameterHandlerBase.cpp.

                                               {
// *************************************
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wuseless-cast"
  // Have the 1 sigma for each parameter in each covariance class, sweet!
  // Don't want to change the prior array because that's what determines our likelihood
  // Want to change the _fPropVal, _fCurrVal, _fPreFitValue
  // _fPreFitValue and the others will already be set
  for (int i = 0; i < _fNumPar; ++i) {
    // Check if parameter is fixed first: if so don't randomly throw
    if (IsParameterFixed(i)) continue;
    // Check that the sigma range is larger than the parameter range
    // If not, throw in the valid parameter range instead
    const double paramrange = _fUpBound[i] - _fLowBound[i];
    const double sigma = sqrt((*covMatrix)(i,i));
    double throwrange = sigma;
    if (paramrange < sigma) throwrange = paramrange;
 
    _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange);
    // Try again if we the initial parameter proposal falls outside of the range of the parameter
    int throws = 0;
    while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
      if (throws > 1000) {
        MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
        MACH3LOG_WARN("Matrix: {}", matrixName);
        MACH3LOG_WARN("Param: {}", _fNames[i]);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange);
      throws++;
    }
    MACH3LOG_INFO("Setting current step in {} param {} = {} from {}", matrixName, i, _fPropVal[i], _fCurrVal[i]);
    _fCurrVal[i] = _fPropVal[i];
  }
  #pragma GCC diagnostic pop
  if (pca) PCAObj->TransferToPCA();
 
  // KS: At the end once we are happy with proposal do special proposal
  SpecialStepProposal();
}

Randomize()

void ParameterHandlerBase::Randomize

(

)

"Randomize" the parameters in the covariance class for the proposed step. Used the proposal kernel and the current parameter value to set proposed step

Definition at line 442 of file ParameterHandlerBase.cpp.

                                                {
// ************************************************
  if (!pca) {
    //KS: By multithreading here we gain at least factor 2 with 8 threads with ND only fit
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < _fNumPar; ++i) {
      // If parameter isn't fixed
      if (!IsParameterFixed(i) > 0.0) {
        randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
        // If parameter IS fixed
      } else {
        randParams[i] = 0.0;
      }
    } // end for
  // If we're in the PCA basis we instead throw parameters there (only _fNumParPCA parameter)
  } else {
    // Scale the random parameters by the sqrt of eigen values for the throw
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < PCAObj->GetNumberPCAedParameters(); ++i)
    {
      // If parameter IS fixed or out of bounds
      if (PCAObj->IsParameterFixedPCA(i)) {
        randParams[i] = 0.0;
      } else {
        randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0,1);
      }
    }
  }
}

ReserveMemory()

void ParameterHandlerBase::ReserveMemory ( const int  size )

protected

Initialise vectors with parameters information.

Parameters

size	integer telling size to which we will resize all vectors/allocate memory

Definition at line 216 of file ParameterHandlerBase.cpp.

                                                          {
// ********************************************
  if (SizeVec <= 0) {
    MACH3LOG_CRITICAL("Covariance matrix {} has {} entries!", GetName(), SizeVec);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  _fNames = std::vector<std::string>(SizeVec);
  _fFancyNames = std::vector<std::string>(SizeVec);
 
  _fPreFitValue     = std::vector<double>(SizeVec, 1.0);
  _fError           = std::vector<double>(SizeVec, 1.0);
  _fCurrVal         = std::vector<double>(SizeVec, 0.0);
  _fPropVal         = std::vector<M3::float_t>(SizeVec, 0.0);
  _fLowBound        = std::vector<double>(SizeVec, -999.99);
  _fUpBound         = std::vector<double>(SizeVec, 999.99);
  _fFlatPrior       = std::vector<bool>(SizeVec, false);
  _fIndivStepScale  = std::vector<double>(SizeVec, 1.0);
 
  corr_throw = new double[SizeVec];
  // set random parameter vector (for correlated steps)
  randParams = new double[SizeVec];
 
  // Set the defaults to true
  for(int i = 0; i < SizeVec; i++) {
    corr_throw[i] = 0.0;
    randParams[i] = 0.0;
  }
 
  InvertCovMatrix.resize(SizeVec, std::vector<double>(SizeVec, 0.0));
  throwMatrixCholDecomp = new double*[SizeVec]();
  // Set the defaults to true
  for(int i = 0; i < SizeVec; i++) {
    throwMatrixCholDecomp[i] = new double[SizeVec]();
    for (int j = 0; j < SizeVec; j++) {
      throwMatrixCholDecomp[i][j] = 0.;
    }
  }
 
  _fGlobalStepScale = 1.0;
}

ResetIndivStepScale()

void ParameterHandlerBase::ResetIndivStepScale

(

)

Adaptive Step Tuning Stuff.

Definition at line 869 of file ParameterHandlerBase.cpp.

                                               {
// ********************************************
  std::vector<double> stepScales(_fNumPar, 1.0);
  _fGlobalStepScale = 1.0;
  SetIndivStepScale(stepScales);
}

RetPointer()

const M3::float_t* ParameterHandlerBase::RetPointer ( const int  iParam ) const

inline

DB Pointer return to param position.

Parameters

iParam

The index of the parameter in the vector.

Returns

A pointer to the parameter value at the specified index.

Warning

ETA - This might be a bit squiffy? If the vector gots moved from say a push_back then the pointer is no longer valid... maybe need a better way to deal with this? It was fine before when the return was to an element of a new array. There must be a clever C++ way to be careful

Definition at line 206 of file ParameterHandlerBase.h.

{return &(_fPropVal.data()[iParam]);}

SanitizeAdaption()

void ParameterHandlerBase::SanitizeAdaption ( ) const

protected

Perform sanity check to ensure adaption isn't misbehaving before fit starts.

Definition at line 966 of file ParameterHandlerBase.cpp.

                                                  {
// ********************************************
  for (size_t i = 0; i < FlipParameterIndex.size(); ++i) {
    const int index = FlipParameterIndex[i];
    if(!param_skip_adapt_flags[index]) {
      MACH3LOG_ERROR("You enabled adaption for parameter which has enabled flipping ({})", _fFancyNames[index]);
      MACH3LOG_ERROR("Right now flipping and adapting doesn't work very well");
      MACH3LOG_ERROR("Please skip adaption for param {}, using ParametersToSkip option in config", _fFancyNames[index]);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  // HH: Loop over correlations to check if any skipped parameter is correlated with adapted one
  // We don't want to change one parameter while keeping the other fixed as this would
  // lead to weird penalty terms in the prior after adapting
  double max_correlation = 0.01; // Define a threshold for significant correlation above which we throw an error
  for (int i = 0; i < _fNumPar; ++i) {
    for (int j = 0; j <= i; ++j) {
      // The symmetry should have been checked during the Init phase
      if(param_skip_adapt_flags[i] && !param_skip_adapt_flags[j]) {
        double corr = (*covMatrix)(i,j)/std::sqrt((*covMatrix)(i,i)*(*covMatrix)(j,j));
        if(std::fabs(corr) > max_correlation) {
          MACH3LOG_ERROR("Correlation between skipped parameter {} ({}) and non-skipped parameter {} ({}) is {:.6e}, above the allowed threshold of {:.6e}.",
                         i, _fFancyNames[i], j, _fFancyNames[j], corr, max_correlation);
          throw MaCh3Exception(__FILE__, __LINE__);
        }
      }
    }
  }
}

SaveAdaptiveToFile()

void ParameterHandlerBase::SaveAdaptiveToFile ( const std::string &  outFileName, const std::string &  systematicName  )

inline

Save adaptive throw matrix to file.

Definition at line 170 of file ParameterHandlerBase.h.

                                                                                         {
    AdaptiveHandler->SaveAdaptiveToFile(outFileName, systematicName); }

SaveUpdatedMatrixConfig()

void ParameterHandlerBase::SaveUpdatedMatrixConfig

(

)

KS: After step scale, prefit etc. value were modified save this modified config.

Definition at line 1196 of file ParameterHandlerBase.cpp.

                                                   {
// ********************************************
  if (!_fYAMLDoc)
  {
    MACH3LOG_CRITICAL("Yaml node hasn't been initialised for matrix {}, something is not right", matrixName);
    MACH3LOG_CRITICAL("I am not throwing error but should be investigated");
    return;
  }
 
  YAML::Node copyNode = _fYAMLDoc;
  int i = 0;
 
  for (YAML::Node param : copyNode["Systematics"])
  {
    //KS: Feel free to update it, if you need updated prefit value etc
    param["Systematic"]["StepScale"]["MCMC"] = M3::Utils::FormatDouble(_fIndivStepScale[i], 4);
    i++;
  }
  // Save the modified node to a file
  std::ofstream fout("Modified_Matrix.yaml");
  fout << copyNode;
  fout.close();
}

SetBranches()

void ParameterHandlerBase::SetBranches	(	TTree &	tree,
		const bool	SaveProposal = false
	)

set branches for output file

Parameters

tree	Tree to which we will save branches
SaveProposal	Normally we only save parameter after is accepted, for debugging purpose it is helpful to see also proposed values. That's what this variable controls

Definition at line 662 of file ParameterHandlerBase.cpp.

                                                                     {
// ********************************************
  // loop over parameters and set a branch
  for (int i = 0; i < _fNumPar; ++i) {
    tree.Branch(_fNames[i].c_str(), &_fCurrVal[i], Form("%s/D", _fNames[i].c_str()));
  }
  // When running PCA, also save PCA parameters
  if (pca) {
    PCAObj->SetBranches(tree, SaveProposal, _fNames);
  }
  if(SaveProposal)
  {
    // loop over parameters and set a branch
    for (int i = 0; i < _fNumPar; ++i) {
      tree.Branch(Form("%s_Prop", _fNames[i].c_str()), &_fPropVal[i], Form("%s_Prop/D", _fNames[i].c_str()));
    }
  }
  if(use_adaptive && AdaptiveHandler->GetUseRobbinsMonro()){
    tree.Branch(Form("GlobalStepScale_%s", GetName().c_str()), &_fGlobalStepScale, Form("GlobalStepScale_%s/D", GetName().c_str()));
  }
}

SetCovMatrix()

void ParameterHandlerBase::SetCovMatrix

(

TMatrixDSym *

cov

)

Set covariance matrix.

Parameters

cov	Covariance matrix which we set and will be used later for evaluation of penalty term

Definition at line 194 of file ParameterHandlerBase.cpp.

                                                        {
// ********************************************
  if (cov == nullptr) {
    MACH3LOG_ERROR("Could not find covariance matrix you provided to {}", __func__ );
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  covMatrix = cov;
 
  invCovMatrix = static_cast<TMatrixDSym *>(cov->Clone());
  invCovMatrix->Invert();
  //KS: ROOT has bad memory management, using standard double means we can decrease most operation by factor 2 simply due to cache hits
  for (int i = 0; i < _fNumPar; i++)
  {
    for (int j = 0; j < _fNumPar; ++j)
    {
      InvertCovMatrix[i][j] = (*invCovMatrix)(i,j);
    }
  }
 
  SetThrowMatrix(cov);
}

SetFixAllParameters()

void ParameterHandlerBase::SetFixAllParameters

(

)

Set all parameters to be fixed at prior values.

Definition at line 716 of file ParameterHandlerBase.cpp.

                                               {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  for (int i = 0; i < _fNumPar; ++i)
    if(!IsParameterFixed(i)) ToggleFixParameter(i);
}

SetFixParameter() [1/2]

void ParameterHandlerBase::SetFixParameter

(

const int

i

)

Set parameter to be fixed at prior value.

Parameters

i	Parameter index

Definition at line 724 of file ParameterHandlerBase.cpp.

                                                      {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  if(!IsParameterFixed(i)) ToggleFixParameter(i);
}

SetFixParameter() [2/2]

void ParameterHandlerBase::SetFixParameter

(

const std::string &

name

)

Set parameter to be fixed at prior value.

Parameters

name	Name of the parameter to be fixed

Definition at line 731 of file ParameterHandlerBase.cpp.

                                                                {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  if(!IsParameterFixed(name)) ToggleFixParameter(name);
}

SetFlatPrior()

void ParameterHandlerBase::SetFlatPrior	(	const int	i,
		const bool	eL
	)

Set if parameter should have flat prior or not.

Parameters

i	Parameter index
eL	bool telling if it will be flat or not

Definition at line 805 of file ParameterHandlerBase.cpp.

                                                                  {
// ********************************************
  if (i > _fNumPar) {
    MACH3LOG_INFO("Can't {} for Cov={}/Param={} because size of Covariance = {}", __func__, GetName(), i, _fNumPar);
    MACH3LOG_ERROR("Fix this in your config file please!");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  } else {
    if(eL){
      MACH3LOG_INFO("Setting {} (parameter {}) to flat prior", GetParName(i), i);
    }
    else{
      // HW :: This is useful
      MACH3LOG_INFO("Setting {} (parameter {}) to non-flat prior", GetParName(i), i);
    }
    _fFlatPrior[i] = eL;
  }
}

SetFreeAllParameters()

void ParameterHandlerBase::SetFreeAllParameters

(

)

Set all parameters to be treated as free.

Definition at line 738 of file ParameterHandlerBase.cpp.

                                                {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  for (int i = 0; i < _fNumPar; ++i)
    if(IsParameterFixed(i)) ToggleFixParameter(i);
}

SetFreeParameter() [1/2]

void ParameterHandlerBase::SetFreeParameter

(

const int

i

)

Set parameter to be treated as free.

Parameters

i	Parameter index

Definition at line 746 of file ParameterHandlerBase.cpp.

                                                       {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  if(IsParameterFixed(i)) ToggleFixParameter(i);
}

SetFreeParameter() [2/2]

void ParameterHandlerBase::SetFreeParameter

(

const std::string &

name

)

Set parameter to be treated as free.

Parameters

name	Name of the parameter to be treated as free

Definition at line 753 of file ParameterHandlerBase.cpp.

                                                                 {
// ********************************************
  // Check if the parameter is fixed and if not, toggle fix it
  if(IsParameterFixed(name)) ToggleFixParameter(name);
}

SetIndivStepScale() [1/2]

void ParameterHandlerBase::SetIndivStepScale ( const int  ParameterIndex, const double  StepScale  )

inline

DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constructors)

Parameters

ParameterIndex	Parameter Index
StepScale	Value of individual step scale

Definition at line 91 of file ParameterHandlerBase.h.

{ _fIndivStepScale.at(ParameterIndex) = StepScale; }

SetIndivStepScale() [2/2]

void ParameterHandlerBase::SetIndivStepScale

(

const std::vector< double > &

stepscale

)

DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constructors)

Parameters

stepscale

Vector of individual step scale, should have same

Definition at line 824 of file ParameterHandlerBase.cpp.

                                                                               {
// ********************************************
  if (static_cast<int>(stepscale.size()) != _fNumPar)
  {
    MACH3LOG_WARN("Stepscale vector not equal to number of parameters. Quitting..");
    MACH3LOG_WARN("Size of argument vector: {}", stepscale.size());
    MACH3LOG_WARN("Expected size: {}", _fNumPar);
    return;
  }
 
  for (int iParam = 0 ; iParam < _fNumPar; iParam++) {
    _fIndivStepScale[iParam] = stepscale[iParam];
  }
  PrintIndivStepScale();
}

SetIndivStepScaleForSkippedAdaptParams()

void ParameterHandlerBase::SetIndivStepScaleForSkippedAdaptParams

(

)

Set individual step scale for parameters which are skipped during adaption to initial values.

Definition at line 877 of file ParameterHandlerBase.cpp.

                                                                  {
// ********************************************
  if (!param_skip_adapt_flags.size()) {
    MACH3LOG_ERROR("Parameter skip adapt flags not set, cannot set individual step scales for skipped parameters.");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  // HH: Cancel the effect of global step scale change for parameters that are not adapting
  for (int i = 0; i <_fNumPar; i++) {
    if (param_skip_adapt_flags[i]) {
      _fIndivStepScale[i] = _fIndivStepScaleInitial[i] * _fGlobalStepScaleInitial / _fGlobalStepScale;
    }
  }
  MACH3LOG_DEBUG("Updating individual step scales for non-adapting parameters to cancel global step scale change.");
  MACH3LOG_DEBUG("Global step scale initial: {}, current: {}", _fGlobalStepScaleInitial, _fGlobalStepScale);
}

SetName()

void ParameterHandlerBase::SetName ( const std::string &  name )

inline

Set matrix name.

Definition at line 41 of file ParameterHandlerBase.h.

{ matrixName = name; }

SetNumberOfSteps()

void ParameterHandlerBase::SetNumberOfSteps ( const int  nsteps )

inline

Set number of MCMC step, when running adaptive MCMC it is updated with given frequency. We need number of steps to determine frequency.

Definition at line 181 of file ParameterHandlerBase.h.

                                          {
    AdaptiveHandler->SetTotalSteps(nsteps);
    if(AdaptiveHandler->AdaptionUpdate()) ResetIndivStepScale();
  }

SetPar()

void ParameterHandlerBase::SetPar	(	const int	i,
		const double	val
	)

Set all the covariance matrix parameters to a user-defined value.

Parameters

i	Parameter index
val	new value which will be set

Definition at line 261 of file ParameterHandlerBase.cpp.

                                                                {
// ********************************************
  MACH3LOG_DEBUG("Over-riding {}: _fPropVal ({}), _fCurrVal ({}), _fPreFitValue ({}) to ({})",
                 GetParFancyName(i), _fPropVal[i], _fCurrVal[i], _fPreFitValue[i], val);
 
  _fPropVal[i] = static_cast<M3::float_t>(val);
  _fCurrVal[i] = val;
  _fPreFitValue[i] = val;
 
  // Transfer the parameter values to the PCA basis
  if (pca) PCAObj->TransferToPCA();
}

SetParameters()

void ParameterHandlerBase::SetParameters

(

const std::vector< double > &

pars = {}

)

Set parameter values using vector, it has to have same size as covariance class.

Parameters

pars	Vector holding new values for every parameter

Definition at line 626 of file ParameterHandlerBase.cpp.

                                                                      {
// ********************************************
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wuseless-cast"
  // If empty, set the proposed to prior
  if (pars.empty()) {
    // For xsec this means setting to the prior (because prior is the prior)
    for (int i = 0; i < _fNumPar; i++) {
      _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i]);
    }
    // If not empty, set the parameters to the specified
  } else {
    if (pars.size() != static_cast<size_t>(_fNumPar)) {
      MACH3LOG_ERROR("Parameter arrays of incompatible size! Not changing parameters! {} has size {} but was expecting {}", matrixName, pars.size(), _fNumPar);
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    int parsSize = int(pars.size());
    for (int i = 0; i < parsSize; i++) {
      //Make sure that you are actually passing a number to set the parameter to
      if(std::isnan(pars[i])) {
        MACH3LOG_ERROR("Trying to set parameter value to a nan for parameter {} in matrix {}. This will not go well!", GetParName(i), matrixName);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      } else {
        _fPropVal[i] = static_cast<M3::float_t>(pars[i]);
      }
    }
  }
  // And if pca make the transfer
  if (pca) {
    PCAObj->TransferToPCA();
    PCAObj->TransferToParam();
  }
  #pragma GCC diagnostic pop
}

SetParCurrProp()

void ParameterHandlerBase::SetParCurrProp	(	const int	i,
		const double	val
	)

Set current parameter value.

Parameters

i	Parameter index
val	new value which will be set

Definition at line 394 of file ParameterHandlerBase.cpp.

                                                                              {
// ********************************************
  _fPropVal[parNo] = static_cast<M3::float_t>(parVal);
  _fCurrVal[parNo] = parVal;
  MACH3LOG_DEBUG("Setting {} (parameter {}) to {})", GetParFancyName(parNo),  parNo, parVal);
  if (pca) PCAObj->TransferToPCA();
}

SetParName()

void ParameterHandlerBase::SetParName ( const int  i, const std::string &  name  )

inline

change parameter name

Parameters

i	Parameter index
name	new name which will be set

Definition at line 45 of file ParameterHandlerBase.h.

{ _fNames.at(i) = name; }

SetParProp()

void ParameterHandlerBase::SetParProp ( const int  i, const double  val  )

inline

Set proposed parameter value.

Parameters

i	Parameter index
val	new value which will be set

Definition at line 59 of file ParameterHandlerBase.h.

                                                 {
    _fPropVal[i] = static_cast<M3::float_t>(val);
    if (pca) PCAObj->TransferToPCA();
  }

SetPrintLength()

void ParameterHandlerBase::SetPrintLength ( const unsigned int  PriLen )

inline

KS: In case someone really want to change this.

Definition at line 96 of file ParameterHandlerBase.h.

{ PrintLength = PriLen; }

SetRandomThrow()

void ParameterHandlerBase::SetRandomThrow ( const int  i, const double  rand  )

inline

Set random value useful for debugging/CI.

Parameters

i	Parameter index
rand	New value for random number

Definition at line 74 of file ParameterHandlerBase.h.

{ randParams[i] = rand;}

SetSingleParameter()

void ParameterHandlerBase::SetSingleParameter	(	const int	parNo,
		const double	parVal
	)

Set value of single param to a given value.

Definition at line 385 of file ParameterHandlerBase.cpp.

                                                                                  {
// *************************************
  _fPropVal[parNo] = static_cast<M3::float_t>(parVal);
  _fCurrVal[parNo] = parVal;
  MACH3LOG_DEBUG("Setting {} (parameter {}) to {})", GetParFancyName(parNo),  parNo, parVal);
  if (pca) PCAObj->TransferToPCA();
}

SetStepScale()

void ParameterHandlerBase::SetStepScale	(	const double	scale,
		const bool	verbose = true
	)

Set global step scale for covariance object.

Parameters

scale	Value of global step scale
verbose	Print that we've changed scale + use warnings [default: true] [25]

Definition at line 685 of file ParameterHandlerBase.cpp.

                                                                              {
// ********************************************
  if(scale <= 0) {
    MACH3LOG_ERROR("You are trying so set StepScale to 0 or negative this will not work");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  if(verbose){
    MACH3LOG_INFO("{} setStepScale() = {}", GetName(), scale);
    const double SuggestedScale = 2.38/std::sqrt(_fNumPar);
    if(std::fabs(scale - SuggestedScale)/SuggestedScale > 1) {
      MACH3LOG_WARN("Defined Global StepScale is {}, while suggested suggested {}", scale, SuggestedScale);
    }
  }
  _fGlobalStepScale = scale;
}

SetSubThrowMatrix()

void ParameterHandlerBase::SetSubThrowMatrix	(	int	first_index,
		int	last_index,
		TMatrixDSym const &	subcov
	)

Definition at line 933 of file ParameterHandlerBase.cpp.

                                                                        {
// ********************************************
  if ((last_index - first_index) >= subcov.GetNrows()) {
    MACH3LOG_ERROR("Trying to SetSubThrowMatrix into range: ({},{}) with a "
                   "submatrix with only {} rows {}",
                   first_index, last_index, subcov.GetNrows(), __func__);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  TMatrixDSym *current_ThrowMatrix =
      static_cast<TMatrixDSym *>(throwMatrix->Clone());
  for (int i = first_index; i <= last_index; ++i) {
    for (int j = first_index; j <= last_index; ++j) {
      current_ThrowMatrix->operator()(i, j) =
          subcov(i - first_index, j - first_index);
    }
  }
 
  SetThrowMatrix(current_ThrowMatrix);
  delete current_ThrowMatrix;
}

SetThrowMatrix()

void ParameterHandlerBase::SetThrowMatrix

(

const TMatrixDSym *

cov

)

Use new throw matrix, used in adaptive MCMC.

Definition at line 895 of file ParameterHandlerBase.cpp.

                                                                {
// ********************************************
   if (cov == nullptr) {
    MACH3LOG_ERROR("Could not find covariance matrix you provided to {}", __func__);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  if (covMatrix->GetNrows() != cov->GetNrows()) {
    MACH3LOG_ERROR("Matrix given for throw Matrix is not the same size as the covariance matrix stored in object!");
    MACH3LOG_ERROR("Stored covariance matrix size: {}", covMatrix->GetNrows());
    MACH3LOG_ERROR("Given matrix size: {}", cov->GetNrows());
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  throwMatrix = static_cast<TMatrixDSym*>(cov->Clone());
  if(use_adaptive && AdaptiveHandler->AdaptionUpdate()) MakeClosestPosDef(throwMatrix);
  else {
    // HW: Prevent spam from adaptive handler
    bool verbose = AdaptiveHandler ? AdaptiveHandler->GetTotalSteps() < 2 : true;
    MakePosDef(throwMatrix, verbose);
  }
 
  auto throwMatrix_CholDecomp = M3::GetCholeskyDecomposedMatrix(*throwMatrix, matrixName);
 
  //KS: ROOT has bad memory management, using standard double means we can decrease most operation by factor 2 simply due to cache hits
  #ifdef MULTITHREAD
  #pragma omp parallel for collapse(2)
  #endif
  for (int i = 0; i < _fNumPar; ++i)
  {
    for (int j = 0; j < _fNumPar; ++j)
    {
      throwMatrixCholDecomp[i][j] = throwMatrix_CholDecomp[i][j];
    }
  }
}

SetThrowMatrixFromFile()

void ParameterHandlerBase::SetThrowMatrixFromFile ( const std::string &  matrix_file_name, const std::string &  matrix_name, const std::string &  means_name  )

protected

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

SetTune()

void ParameterHandlerBase::SetTune

(

const std::string &

TuneName

)

KS: Set proposed parameter values vector to be base on tune values, for example set proposed values to be of generated or maybe PostND.

Definition at line 1222 of file ParameterHandlerBase.cpp.

                                                            {
// ********************************************
  if(Tunes == nullptr) {
    MACH3LOG_ERROR("Tunes haven't been initialised, which are being loaded from YAML, have you used some deprecated constructor");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  auto Values = Tunes->GetTune(TuneName);
 
  SetParameters(Values);
}

SpecialStepProposal()

void ParameterHandlerBase::SpecialStepProposal ( )

protected

Perform Special Step Proposal.

Warning

KS: Following Asher comment we do "Step->Circular Bounds->Flip"

Warning

KS: Following Asher comment we do "Step->Circular Bounds->Flip"

Definition at line 419 of file ParameterHandlerBase.cpp.

                                               {
// ************************************************
 
  // HW It should now automatically set dcp to be with [-pi, pi]
  for (size_t i = 0; i < CircularBoundsIndex.size(); ++i) {
    const int index = CircularBoundsIndex[i];
    if(!IsParameterFixed(index))
      CircularParBounds(index, CircularBoundsValues[i].first, CircularBoundsValues[i].second);
  }
 
  // Okay now we've done the standard steps, we can add in our nice flips hierarchy flip first
  for (size_t i = 0; i < FlipParameterIndex.size(); ++i) {
    const int index = FlipParameterIndex[i];
    if(!IsParameterFixed(index))
      FlipParameterValue(FlipParameterIndex[i], FlipParameterPoint[i]);
  }
}

ThrowParameters()

void ParameterHandlerBase::ThrowParameters

(

)

Throw the parameters according to the covariance matrix. This shouldn't be used in MCMC code ase it can break Detailed Balance;.

Definition at line 285 of file ParameterHandlerBase.cpp.

                                           {
// *************************************
  // First draw new randParams
  Randomize();
 
  M3::MatrixVectorMulti(corr_throw, throwMatrixCholDecomp, randParams, _fNumPar);
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wuseless-cast"
  // KS: We use PCA very rarely on top PCA functionality isn't implemented for this function.
  // Use __builtin_expect to give compiler a hint which option is more likely, which should help
  // with better optimisation. This isn't critical but more to have example
  if (__builtin_expect(!pca, 1)) {
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < _fNumPar; ++i) {
      // Check if parameter is fixed first: if so don't randomly throw
      if (IsParameterFixed(i)) continue;
      _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i] + corr_throw[i]);
 
      int throws = 0;
      // Try again if we the initial parameter proposal falls outside of the range of the parameter
      while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
        randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
        const double corr_throw_single = M3::MatrixVectorMultiSingle(throwMatrixCholDecomp, randParams, _fNumPar, i);
        _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i] + corr_throw_single);
        if (throws > 10000)
        {
          //KS: Since we are multithreading there is danger that those messages
          //will be all over the place, small price to pay for faster code
          MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
          MACH3LOG_WARN("Matrix: {}", matrixName);
          MACH3LOG_WARN("Param: {}", _fNames[i]);
          MACH3LOG_WARN("Setting _fPropVal:  {} to {}", _fPropVal[i], _fPreFitValue[i]);
          MACH3LOG_WARN("I live at {}:{}", __FILE__, __LINE__);
          _fPropVal[i] = static_cast<M3::float_t>(_fPreFitValue[i]);
          //throw MaCh3Exception(__FILE__ , __LINE__ );
        }
        throws++;
      }
      _fCurrVal[i] = _fPropVal[i];
    }
  }
  else
  {
    PCAObj->ThrowParameters(random_number, throwMatrixCholDecomp,
                            randParams, corr_throw,
                            _fPreFitValue, _fLowBound, _fUpBound, _fNumPar);
  } // end if pca
  #pragma GCC diagnostic pop
  // KS: At the end once we are happy with proposal do special proposal
  SpecialStepProposal();
}

ToggleFixParameter() [1/2]

void ParameterHandlerBase::ToggleFixParameter ( const int  i )

private

Toggle fixing parameter at prior values.

Parameters

i	Parameter index

Definition at line 760 of file ParameterHandlerBase.cpp.

                                                         {
// ********************************************
  if(!pca) {
    if (i > _fNumPar) {
      MACH3LOG_ERROR("Can't {} for parameter {} because size of covariance ={}", __func__, i, _fNumPar);
      MACH3LOG_ERROR("Fix this in your config file please!");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    } else {
      _fError[i] *= -1.0;
      if(IsParameterFixed(i)) MACH3LOG_INFO("Setting {}(parameter {}) to fixed at {}", GetParFancyName(i), i, _fCurrVal[i]);
      else MACH3LOG_INFO("Setting {}(parameter {}) free", GetParFancyName(i), i);
    }
    if( (_fCurrVal[i] > _fUpBound[i] || _fCurrVal[i] < _fLowBound[i]) && IsParameterFixed(i) ) {
      MACH3LOG_ERROR("Parameter {} (index {}) is fixed at {}, which is outside of its bounds [{}, {}]", GetParFancyName(i), i, _fCurrVal[i], _fLowBound[i], _fUpBound[i]);
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  } else {
    PCAObj->ToggleFixParameter(i, _fNames);
  }
}

ToggleFixParameter() [2/2]

void ParameterHandlerBase::ToggleFixParameter ( const std::string &  name )

private

Toggle fixing parameter at prior values.

Parameters

name	Name of parameter you want to fix

Definition at line 782 of file ParameterHandlerBase.cpp.

                                                                   {
// ********************************************
  const int Index = GetParIndex(name);
  if(Index != M3::_BAD_INT_) {
    ToggleFixParameter(Index);
    return;
  }
  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore will not fix it", name);
}

UpdateAdaptiveCovariance()

void ParameterHandlerBase::UpdateAdaptiveCovariance

(

)

Method to update adaptive MCMC [15].

Need to adjust the scale every step

Definition at line 1075 of file ParameterHandlerBase.cpp.

                                                    {
// ********************************************
  // Updates adaptive matrix
  // First we update the total means
 
  // Skip this if we're at a large number of steps
  if(AdaptiveHandler->SkipAdaption()) {
    AdaptiveHandler->IncrementNSteps();
    return;
  }
 
  if(AdaptiveHandler->GetUseRobbinsMonro()){
    bool verbose=false;
    #ifdef MACH3_DEBUG
    verbose=true;
    #endif
    AdaptiveHandler->UpdateRobbinsMonroScale();
    SetStepScale(AdaptiveHandler->GetAdaptionScale(), verbose);
    SetIndivStepScaleForSkippedAdaptParams();
  }
 
  // Call main adaption function
  AdaptiveHandler->UpdateAdaptiveCovariance();
 
  // Set scales to 1 * optimal scale
  if(AdaptiveHandler->IndivStepScaleAdapt()) {
    ResetIndivStepScale();
    SetStepScale(AdaptiveHandler->GetAdaptionScale());
    SetIndivStepScaleForSkippedAdaptParams();
  }
 
  if(AdaptiveHandler->UpdateMatrixAdapt()) {
    TMatrixDSym* update_matrix = static_cast<TMatrixDSym*>(AdaptiveHandler->GetAdaptiveCovariance()->Clone());
    UpdateThrowMatrix(update_matrix); //Now we update and continue!
    //Also Save the adaptive to file
    AdaptiveHandler->SaveAdaptiveToFile(AdaptiveHandler->GetOutFileName(), GetName());
  }
 
  AdaptiveHandler->IncrementNSteps();
}

UpdateThrowMatrix()

void ParameterHandlerBase::UpdateThrowMatrix

(

TMatrixDSym *

cov

)

Replaces old throw matrix with new one.

Definition at line 957 of file ParameterHandlerBase.cpp.

                                                             {
// ********************************************
  delete throwMatrix;
  throwMatrix = nullptr;
  SetThrowMatrix(cov);
}

Member Data Documentation

_fCurrVal

std::vector<double> ParameterHandlerBase::_fCurrVal

protected

Current value of the parameter.

Definition at line 424 of file ParameterHandlerBase.h.

_fError

std::vector<double> ParameterHandlerBase::_fError

protected

Prior error on the parameter.

Definition at line 428 of file ParameterHandlerBase.h.

_fFancyNames

std::vector<std::string> ParameterHandlerBase::_fFancyNames

protected

Fancy name for example rather than param_0 it is MAQE, useful for human reading.

Definition at line 416 of file ParameterHandlerBase.h.

_fFlatPrior

std::vector<bool> ParameterHandlerBase::_fFlatPrior

protected

Whether to apply flat prior or not.

Definition at line 436 of file ParameterHandlerBase.h.

_fGlobalStepScale

double ParameterHandlerBase::_fGlobalStepScale

protected

Global step scale applied to all params in this class.

Definition at line 408 of file ParameterHandlerBase.h.

_fGlobalStepScaleInitial

double ParameterHandlerBase::_fGlobalStepScaleInitial

protected

Backup of _fGlobalStepScale for parameters which are skipped during adaption.

Definition at line 442 of file ParameterHandlerBase.h.

_fIndivStepScale

std::vector<double> ParameterHandlerBase::_fIndivStepScale

protected

Individual step scale used by MCMC algorithm.

Definition at line 434 of file ParameterHandlerBase.h.

_fIndivStepScaleInitial

std::vector<double> ParameterHandlerBase::_fIndivStepScaleInitial

protected

Backup of _fIndivStepScale for parameters which are skipped during adaption.

Definition at line 439 of file ParameterHandlerBase.h.

_fLowBound

std::vector<double> ParameterHandlerBase::_fLowBound

protected

Lowest physical bound, parameter will not be able to go beyond it.

Definition at line 430 of file ParameterHandlerBase.h.

_fNames

std::vector<std::string> ParameterHandlerBase::_fNames

protected

ETA _fNames is set automatically in the covariance class to be something like param_i, this is currently to make things compatible with the Diagnostic tools.

Definition at line 414 of file ParameterHandlerBase.h.

_fNumPar

int ParameterHandlerBase::_fNumPar

protected

Number of systematic parameters.

Definition at line 420 of file ParameterHandlerBase.h.

_fPreFitValue

std::vector<double> ParameterHandlerBase::_fPreFitValue

protected

Parameter value dictated by the prior model. Based on it penalty term is calculated.

Definition at line 422 of file ParameterHandlerBase.h.

_fPropVal

std::vector<M3::float_t> ParameterHandlerBase::_fPropVal

protected

Proposed value of the parameter.

Definition at line 426 of file ParameterHandlerBase.h.

_fUpBound

std::vector<double> ParameterHandlerBase::_fUpBound

protected

Upper physical bound, parameter will not be able to go beyond it.

Definition at line 432 of file ParameterHandlerBase.h.

_fYAMLDoc

YAML::Node ParameterHandlerBase::_fYAMLDoc

protected

Stores config describing systematics.

Definition at line 418 of file ParameterHandlerBase.h.

AdaptiveHandler

std::unique_ptr<AdaptiveMCMCHandler> ParameterHandlerBase::AdaptiveHandler

protected

Struct containing information about adaption.

Definition at line 460 of file ParameterHandlerBase.h.

CircularBoundsIndex

std::vector<int> ParameterHandlerBase::CircularBoundsIndex

protected

Indices of parameters with circular bounds.

Definition at line 469 of file ParameterHandlerBase.h.

CircularBoundsValues

std::vector<std::pair<double,double> > ParameterHandlerBase::CircularBoundsValues

protected

Circular bounds for each parameter (lower, upper)

Definition at line 471 of file ParameterHandlerBase.h.

corr_throw

double* ParameterHandlerBase::corr_throw

protected

Result of multiplication of Cholesky matrix and randParams.

Definition at line 406 of file ParameterHandlerBase.h.

covMatrix

TMatrixDSym* ParameterHandlerBase::covMatrix

protected

The covariance matrix.

Definition at line 394 of file ParameterHandlerBase.h.

doSpecialStepProposal

bool ParameterHandlerBase::doSpecialStepProposal

protected

Check if any of special step proposal were enabled.

Definition at line 386 of file ParameterHandlerBase.h.

FlipParameterIndex

std::vector<int> ParameterHandlerBase::FlipParameterIndex

protected

Indices of parameters with flip symmetry.

Definition at line 465 of file ParameterHandlerBase.h.

FlipParameterPoint

std::vector<double> ParameterHandlerBase::FlipParameterPoint

protected

Central points around which parameters are flipped.

Definition at line 467 of file ParameterHandlerBase.h.

inputFile

std::string ParameterHandlerBase::inputFile

protected

The input root file we read in.

Definition at line 389 of file ParameterHandlerBase.h.

invCovMatrix

TMatrixDSym* ParameterHandlerBase::invCovMatrix

protected

The inverse covariance matrix.

Definition at line 396 of file ParameterHandlerBase.h.

InvertCovMatrix

std::vector<std::vector<double> > ParameterHandlerBase::InvertCovMatrix

protected

KS: Same as above but much faster as TMatrixDSym cache miss.

Definition at line 398 of file ParameterHandlerBase.h.

matrixName

std::string ParameterHandlerBase::matrixName

protected

Name of cov matrix.

Definition at line 392 of file ParameterHandlerBase.h.

param_skip_adapt_flags

std::vector<bool> ParameterHandlerBase::param_skip_adapt_flags

protected

Flags telling if parameter should be skipped during adaption.

Definition at line 445 of file ParameterHandlerBase.h.

pca

bool ParameterHandlerBase::pca

protected

perform PCA or not

Definition at line 453 of file ParameterHandlerBase.h.

PCAObj

std::unique_ptr<PCAHandler> ParameterHandlerBase::PCAObj

protected

Struct containing information about PCA.

Definition at line 458 of file ParameterHandlerBase.h.

PrintLength

int ParameterHandlerBase::PrintLength

protected

KS: This is used when printing parameters, sometimes we have super long parameters name, we want to flexibly adjust couts.

Definition at line 411 of file ParameterHandlerBase.h.

random_number

std::vector<std::unique_ptr<TRandom3> > ParameterHandlerBase::random_number

protected

KS: Set Random numbers for each thread so each thread has different seed.

Definition at line 401 of file ParameterHandlerBase.h.

randParams

double* ParameterHandlerBase::randParams

protected

Random number taken from gaussian around prior error used for corr_throw.

Definition at line 404 of file ParameterHandlerBase.h.

throwMatrix

TMatrixDSym* ParameterHandlerBase::throwMatrix

protected

Matrix which we use for step proposal before Cholesky decomposition (not actually used for step proposal)

Definition at line 448 of file ParameterHandlerBase.h.

throwMatrixCholDecomp

double** ParameterHandlerBase::throwMatrixCholDecomp

protected

Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss.

Definition at line 450 of file ParameterHandlerBase.h.

Tunes

std::unique_ptr<ParameterTunes> ParameterHandlerBase::Tunes

protected

Struct containing information about adaption.

Definition at line 462 of file ParameterHandlerBase.h.

use_adaptive

bool ParameterHandlerBase::use_adaptive

protected

Are we using AMCMC?

Definition at line 455 of file ParameterHandlerBase.h.

---

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

• Parameters/ParameterHandlerBase.h
• Parameters/ParameterHandlerBase.cpp