ParameterHandlerBase Class Reference

Base class responsible for handling of systematic error parameters. Capable of using PCA or using adaptive throw matrix. More...

#include <Parameters/ParameterHandlerBase.h>

Inheritance diagram for ParameterHandlerBase:

Collaboration diagram for ParameterHandlerBase:

Public Member Functions
	ParameterHandlerBase (const std::vector< std::string > &YAMLFile, std::string name, double threshold=-1, int FirstPCAdpar=-999, int LastPCAdpar=-999)
	ETA - constructor for a YAML file. More...
	ParameterHandlerBase (std::string name, std::string file, double threshold=-1, int FirstPCAdpar=-999, int LastPCAdpar=-999)
	"Usual" constructors from root file More...
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...
std::string	GetInputFile () const
	Get name of input file. 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 (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 ()
	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)
	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	GetParInit (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	PrintNominal () const
	Print prior value for every parameter. More...
void	PrintNominalCurrProp () const
	Print prior, current and proposed value for each parameter. More...
void	PrintParameters () const
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...
void	ToggleFixAllParameters ()
	Toggle fixing parameters at prior values. More...
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...
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	Init (const std::string &name, const std::string &file)
	Initialisation of the class using matrix from root file. More...
void	Init (const std::vector< std::string > &YAMLFile)
	Initialisation of the class using config. More...
void	ReserveMemory (const int size)
	Initialise vectors with parameters information. More...
void	MakePosDef (TMatrixDSym *cov=nullptr)
	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...
bool	AppliesToSample (const int SystIndex, const std::string &SampleName) const
	Check if parameter is affecting given sample name. 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...
const 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< std::vector< std::string > >	_fSampleNames
	Tells to which samples object param should be applied. 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...

Detailed Description

Base class responsible for handling of systematic error parameters. Capable of using PCA or using adaptive throw matrix.

Author

Dan Barrow

Ed Atkin

Kamil Skwarczynski

Definition at line 15 of file ParameterHandlerBase.h.

Constructor & Destructor Documentation

ParameterHandlerBase() [1/2]

ParameterHandlerBase::ParameterHandlerBase	(	const std::vector< std::string > &	YAMLFile,
		std::string	name,
		double	threshold = -1,
		int	FirstPCAdpar = -999,
		int	LastPCAdpar = -999
	)

ETA - constructor for a YAML file.

Parameters

YAMLFile	A vector of strings representing the YAML files used for initialisation of matrix
name	Matrix name
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.

Definition at line 27 of file ParameterHandlerBase.cpp.

                     : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true) {
// ********************************************
  MACH3LOG_INFO("Constructing instance of ParameterHandler using");
  doSpecialStepProposal = false;
  for(unsigned int i = 0; i < YAMLFile.size(); i++)
  {
    MACH3LOG_INFO("{}", YAMLFile[i]);
  }
  MACH3LOG_INFO("as an input");
 
  if (threshold < 0 || threshold >= 1) {
    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;
  }
 
  Init(YAMLFile);
  // Call the innocent helper function
  if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
}

ParameterHandlerBase() [2/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 8 of file ParameterHandlerBase.cpp.

                     : inputFile(file), pca(true) {
// ********************************************
  MACH3LOG_DEBUG("Constructing instance of ParameterHandler");
  doSpecialStepProposal = 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;
  }
  Init(name, file);
 
  // Call the innocent helper function
  if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
}

~ParameterHandlerBase()

ParameterHandlerBase::~ParameterHandlerBase ( )

virtual

Destructor.

Definition at line 53 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 754 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();
  }
}

AppliesToSample()

bool ParameterHandlerBase::AppliesToSample ( const int  SystIndex, const std::string &  SampleName  ) const

protected

Check if parameter is affecting given sample name.

Parameters

SystIndex	number of parameter
SampleName	The Sample name used to filter parameters.

Definition at line 1468 of file ParameterHandlerBase.cpp.

                                                                                                 {
// ********************************************
  // Empty means apply to all
  if (_fSampleNames[SystIndex].size() == 0) return true;
 
  // Make a copy and to lower case to not be case sensitive
  std::string SampleNameCopy = SampleName;
  std::transform(SampleNameCopy.begin(), SampleNameCopy.end(), SampleNameCopy.begin(), ::tolower);
 
  // Check for unsupported wildcards in SampleNameCopy
  if (SampleNameCopy.find('*') != std::string::npos) {
    MACH3LOG_ERROR("Wildcards ('*') are not supported in sample name: '{}'", SampleName);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  bool Applies = false;
 
  for (size_t i = 0; i < _fSampleNames[SystIndex].size(); i++) {
    // Convert to low case to not be case sensitive
    std::string pattern = _fSampleNames[SystIndex][i];
    std::transform(pattern.begin(), pattern.end(), pattern.begin(), ::tolower);
 
    // Replace '*' in the pattern with '.*' for regex matching
    std::string regexPattern = "^" + std::regex_replace(pattern, std::regex("\\*"), ".*") + "$";
    try {
      std::regex regex(regexPattern);
      if (std::regex_match(SampleNameCopy, regex)) {
        Applies = true;
        break;
      }
    } catch (const std::regex_error& e) {
      // Handle regex error (for invalid patterns)
      MACH3LOG_ERROR("Regex error: {}", e.what());
    }
  }
  return Applies;
}

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 827 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 855 of file ParameterHandlerBase.cpp.

                                                       {
// ********************************************
  int NOutside = 0;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:NOutside)
  #endif
  for (int i = 0; i < _fNumPar; ++i){
    if(_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]){
      NOutside++;
    }
  }
  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 777 of file ParameterHandlerBase.cpp.

                                                                                                         {
// *************************************
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wuseless-cast"
  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));
  }
  #pragma GCC diagnostic pop
}

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 68 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 729 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 381 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__);
    }
  }
 
  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 790 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 320 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 317 of file ParameterHandlerBase.h.

{ return _fYAMLDoc; }

GetCorrelationMatrix()

TH2D * ParameterHandlerBase::GetCorrelationMatrix

(

)

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 1415 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 131 of file ParameterHandlerBase.h.

{ return corr_throw[i]; }

GetCovMatrix()

TMatrixDSym* ParameterHandlerBase::GetCovMatrix ( ) const

inline

Return covariance matrix.

Definition at line 123 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 154 of file ParameterHandlerBase.h.

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

GetDoAdaption()

bool ParameterHandlerBase::GetDoAdaption ( ) const

inline

Do we adapt or not.

Definition at line 173 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 157 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 135 of file ParameterHandlerBase.h.

{ return _fFlatPrior[i]; }

GetGlobalStepScale()

double ParameterHandlerBase::GetGlobalStepScale ( ) const

inline

Get global step scale for covariance object.

Definition at line 239 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 237 of file ParameterHandlerBase.h.

{return _fIndivStepScale.at(ParameterIndex); }

GetInputFile()

std::string ParameterHandlerBase::GetInputFile ( ) const

inline

Get name of input file.

Definition at line 150 of file ParameterHandlerBase.h.

{ return inputFile; }

GetInvCovMatrix() [1/2]

TMatrixDSym* ParameterHandlerBase::GetInvCovMatrix ( ) const

inline

Return inverted covariance matrix.

Definition at line 125 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 127 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 870 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 234 of file ParameterHandlerBase.h.

{ return _fLowBound[i]; }

GetName()

std::string ParameterHandlerBase::GetName ( ) const

inline

Get name of covariance.

Definition at line 138 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 242 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 212 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 222 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 224 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 148 of file ParameterHandlerBase.h.

{return _fFancyNames[i];}

GetParIndex()

int ParameterHandlerBase::GetParIndex

(

const std::string &

name

)

const

Get index based on name.

Definition at line 959 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;
}

GetParInit()

double ParameterHandlerBase::GetParInit ( const int  i ) const

inline

Get prior parameter value.

Parameters

i	Parameter index

Definition at line 228 of file ParameterHandlerBase.h.

{ return _fPreFitValue[i]; }

GetParName()

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

inline

Get name of parameter.

Parameters

i	Parameter index

Definition at line 141 of file ParameterHandlerBase.h.

{return _fNames[i];}

GetParProp()

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

inline

Get proposed parameter value.

Parameters

i	Parameter index

Definition at line 219 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 209 of file ParameterHandlerBase.h.

{return _fPropVal;}

GetPCAHandler()

PCAHandler* ParameterHandlerBase::GetPCAHandler ( ) const

inline

Get pointer for PCAHandler.

Definition at line 332 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 214 of file ParameterHandlerBase.h.

{return _fPreFitValue;}

GetProposed()

std::vector< double > ParameterHandlerBase::GetProposed

(

)

const

Get vector of all proposed parameter values.

Definition at line 526 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 74 of file ParameterHandlerBase.h.

{ return randParams[i];}

GetThrowMatrix() [1/2]

TMatrixDSym* ParameterHandlerBase::GetThrowMatrix ( ) const

inline

Get matrix used for step proposal.

Definition at line 186 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 188 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 231 of file ParameterHandlerBase.h.

{ return _fUpBound[i];}

Init() [1/2]

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

protected

Initialisation of the class using matrix from root file.

Definition at line 95 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));
  }
  // Not using adaptive by default
  use_adaptive = false;
  // Set the covariance matrix
  _fNumPar = CovMat->GetNrows();
 
  InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
  throwMatrixCholDecomp = new double*[_fNumPar]();
  // Set the defaults to true
  for(int i = 0; i < _fNumPar; i++) {
    throwMatrixCholDecomp[i] = new double[_fNumPar]();
    for (int j = 0; j < _fNumPar; j++) {
      throwMatrixCholDecomp[i][j] = 0.;
    }
  }
  SetName(name);
  MakePosDef(CovMat);
  SetCovMatrix(CovMat);
  if (_fNumPar <= 0) {
    MACH3LOG_CRITICAL("Covariance matrix {} has {} entries!", GetName(), _fNumPar);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  ReserveMemory(_fNumPar);
 
  infile->Close();
 
  MACH3LOG_INFO("Created covariance matrix named: {}", GetName());
  MACH3LOG_INFO("from file: {}", file);
  delete infile;
}

Init() [2/2]

void ParameterHandlerBase::Init ( const std::vector< std::string > &  YAMLFile )

protected

Initialisation of the class using config.

Parameters

YAMLFile

A vector of strings representing the YAML files used for initialisation of matrix

Definition at line 156 of file ParameterHandlerBase.cpp.

                                                                    {
// ********************************************
  std::map<std::pair<int, int>, std::unique_ptr<TMatrixDSym>> ThrowSubMatrixOverrides;
  int running_num_file_pars = 0;
 
  _fYAMLDoc["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
  for(unsigned int i = 0; i < YAMLFile.size(); i++)
  {
    YAML::Node YAMLDocTemp = M3OpenConfig(YAMLFile[i]);
 
    if (YAMLDocTemp["ThrowMatrixOverride"]) { // LP: this allows us to put in
                                              // proposal matrix overrides per
                                              // parameter-containing file, add
                                              // the block diagonal proposal
                                              // matrix to a list and overwrite
                                              // the throw matrix after set up.
      auto filename =
          YAMLDocTemp["ThrowMatrixOverride"]["file"].as<std::string>();
      TFile *submatrix_file = TFile::Open(filename.c_str());
 
      auto matrixname =
          YAMLDocTemp["ThrowMatrixOverride"]["matrix"].as<std::string>();
      std::unique_ptr<TMatrixDSym> submatrix{
          submatrix_file->Get<TMatrixDSym>(matrixname.c_str())};
      if (!submatrix) {
        MACH3LOG_CRITICAL("Covariance matrix {} doesn't exist in file: {}",
                          matrixname, filename);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      auto numrows = submatrix->GetNrows();
      // LP: the -1 here is because we specify the last index for consistency
      // with PCAHandler, not the first index after the end as is more common
      // throughout computer science...
      ThrowSubMatrixOverrides[{running_num_file_pars,
                               running_num_file_pars + (numrows - 1)}] =
          std::move(submatrix);
 
      // LP: check names by default, but have option to disable check if you
      // know what you're doing
      if (!bool(YAMLDocTemp["ThrowMatrixOverride"]["check_names"]) ||
          YAMLDocTemp["ThrowMatrixOverride"]["check_names"].as<bool>()) {
        auto nametree = submatrix_file->Get<TTree>("param_names");
        if (!nametree) {
          MACH3LOG_CRITICAL("TTree param_names doesn't exist in file: {}. Set "
                            "ThrowMatrixOverride: {{ check_names: False }} to "
                            "disable this check.",
                            filename);
          throw MaCh3Exception(__FILE__, __LINE__);
        }
        std::string *param_name = nullptr;
        nametree->SetBranchAddress("name", &param_name);
 
        if (nametree->GetEntries() != int(YAMLDocTemp["Systematics"].size())) {
          MACH3LOG_CRITICAL("TTree param_names in file: {} has {} entries, but "
                            "the corresponding yaml file only declares {} "
                            "parameters. Set ThrowMatrixOverride: {{ "
                            "check_names: False }} to disable this check.",
                            filename, nametree->GetEntries(),
                            YAMLDocTemp["Systematics"].size());
          throw MaCh3Exception(__FILE__, __LINE__);
        }
 
        int pit = 0;
        for (const auto &param : YAMLDocTemp["Systematics"]) {
          nametree->GetEntry(pit++);
          auto yaml_pname = Get<std::string>(
              param["Systematic"]["Names"]["FancyName"], __FILE__, __LINE__);
          if ((*param_name) != yaml_pname) {
            MACH3LOG_CRITICAL(
                "TTree param_names in file: {} at entry {} has parameter {}, "
                "but "
                "the corresponding yaml parameter is named {}. Set "
                "ThrowMatrixOverride: {{ "
                "check_names: False }} to disable this check.",
                filename, pit, (*param_name), yaml_pname);
            throw MaCh3Exception(__FILE__, __LINE__);
          }
        }
      }
      submatrix_file->Close();
    }
 
    for (const auto& item : YAMLDocTemp["Systematics"]) {
      _fYAMLDoc["Systematics"].push_back(item);
      running_num_file_pars++;
    }
  }
 
  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));
  }
  PrintLength = 35;
 
  // Set the covariance matrix
  _fNumPar = int(_fYAMLDoc["Systematics"].size());
 
  use_adaptive = false;
 
  InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
  throwMatrixCholDecomp = new double*[_fNumPar]();
  for(int i = 0; i < _fNumPar; i++) {
    throwMatrixCholDecomp[i] = new double[_fNumPar]();
    for (int j = 0; j < _fNumPar; j++) {
      throwMatrixCholDecomp[i][j] = 0.;
    }
  }
  ReserveMemory(_fNumPar);
 
  TMatrixDSym* _fCovMatrix = new TMatrixDSym(_fNumPar);
  int i = 0;
  std::vector<std::map<std::string,double>> Correlations(_fNumPar);
  std::map<std::string, int> CorrNamesMap;
 
  //ETA - read in the systematics. Would be good to add in some checks to make sure
  //that there are the correct number of entries i.e. are the _fNumPar for Names,
  //PreFitValues etc etc.
 
  for (auto const &param : _fYAMLDoc["Systematics"])
  {
    _fFancyNames[i] = Get<std::string>(param["Systematic"]["Names"]["FancyName"], __FILE__ , __LINE__);
    _fPreFitValue[i] = Get<double>(param["Systematic"]["ParameterValues"]["PreFitValue"], __FILE__ , __LINE__);
    _fIndivStepScale[i] = Get<double>(param["Systematic"]["StepScale"]["MCMC"], __FILE__ , __LINE__);
    _fError[i] = Get<double>(param["Systematic"]["Error"], __FILE__ , __LINE__);
    _fSampleNames[i] = GetFromManager<std::vector<std::string>>(param["Systematic"]["SampleNames"], {}, __FILE__, __LINE__);
    if(_fError[i] <= 0) {
      MACH3LOG_ERROR("Error for param {}({}) is negative and equal to {}", _fFancyNames[i], i, _fError[i]);
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    //ETA - a bit of a fudge but works
    auto TempBoundsVec = GetBounds(param["Systematic"]["ParameterBounds"]);
    _fLowBound[i] = TempBoundsVec[0];
    _fUpBound[i] = TempBoundsVec[1];
 
    //ETA - now for parameters which are optional and have default values
    _fFlatPrior[i] = GetFromManager<bool>(param["Systematic"]["FlatPrior"], false, __FILE__ , __LINE__);
 
    // Allow to fix param, this setting should be used only for params which are permanently fixed like baseline, please use global config for fixing param more flexibly
    if(GetFromManager<bool>(param["Systematic"]["FixParam"], false, __FILE__ , __LINE__)) {
      ToggleFixParameter(_fFancyNames[i]);
    }
 
    if(param["Systematic"]["SpecialProposal"]) {
      EnableSpecialProposal(param["Systematic"]["SpecialProposal"], i);
    }
 
    //Fill the map to get the correlations later as well
    CorrNamesMap[param["Systematic"]["Names"]["FancyName"].as<std::string>()]=i;
 
    //Also loop through the correlations
    if(param["Systematic"]["Correlations"]) {
      for(unsigned int Corr_i = 0; Corr_i < param["Systematic"]["Correlations"].size(); ++Corr_i){
        for (YAML::const_iterator it = param["Systematic"]["Correlations"][Corr_i].begin(); it!=param["Systematic"]["Correlations"][Corr_i].end();++it) {
          Correlations[i][it->first.as<std::string>()] = it->second.as<double>();
        }
      }
    }
    i++;
  } // end loop over para
  if(i != _fNumPar) {
    MACH3LOG_CRITICAL("Inconsistent number of params in Yaml  {} vs {}, this indicate wrong syntax", i, i, _fNumPar);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  // ETA Now that we've been through all systematic let's fill the covmatrix
  //This makes the root TCov from YAML
  for(int j = 0; j < _fNumPar; j++) {
    (*_fCovMatrix)(j, j) = _fError[j]*_fError[j];
    //Get the map of parameter name to correlation from the Correlations object
    for (auto const& pair : Correlations[j]) {
      auto const& key = pair.first;
      auto const& val = pair.second;
      int index = -1;
      //If you found the parameter name then get the index
      if (CorrNamesMap.find(key) != CorrNamesMap.end()) {
        index = CorrNamesMap[key];
      } else {
        MACH3LOG_ERROR("Parameter {} not in list! Check your spelling?", key);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      double Corr1 = val;
      double Corr2 = 0;
      if(Correlations[index].find(_fFancyNames[j]) != Correlations[index].end()) {
        Corr2 = Correlations[index][_fFancyNames[j]];
        //Do they agree to better than float precision?
        if(std::abs(Corr2 - Corr1) > FLT_EPSILON) {
          MACH3LOG_ERROR("Correlations are not equal between {} and {}", _fFancyNames[j], key);
          MACH3LOG_ERROR("Got : {} and {}", Corr2, Corr1);
          throw MaCh3Exception(__FILE__ , __LINE__ );
        }
      } else {
        MACH3LOG_ERROR("Correlation does not appear reciprocally between {} and {}", _fFancyNames[j], key);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      (*_fCovMatrix)(j, index)= (*_fCovMatrix)(index, j) = Corr1*_fError[j]*_fError[index];
    }
  }
 
  //Now make positive definite
  MakePosDef(_fCovMatrix);
  SetCovMatrix(_fCovMatrix);
 
  if (_fNumPar <= 0) {
    MACH3LOG_ERROR("ParameterHandler object has {} systematics!", _fNumPar);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  for(auto const & matovr : ThrowSubMatrixOverrides){
    SetSubThrowMatrix(matovr.first.first, matovr.first.second, *matovr.second);
  }
 
  Tunes = std::make_unique<ParameterTunes>(_fYAMLDoc["Systematics"]);
 
  MACH3LOG_INFO("Created covariance matrix from files: ");
  for(const auto &file : YAMLFile){
    MACH3LOG_INFO("{} ", file);
  }
  MACH3LOG_INFO("----------------");
  MACH3LOG_INFO("Found {} systematics parameters in total", _fNumPar);
  MACH3LOG_INFO("----------------");
}

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 1226 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
  std::vector<std::string> params_to_skip = GetFromManager<std::vector<std::string>>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ParametersToSkip"], {});
  // Build a list of skip flags
  param_skip_adapt_flags.resize(_fNumPar, false);
  for (int i = 0; i <_fNumPar; ++i) {
    for (const auto& name : params_to_skip) {
      if(name == _fFancyNames[i]) {
        param_skip_adapt_flags[i] = true;
        break;
      }
    }
  }
  // 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__);
        }
      }
    }
  }
  // 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) {
    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 298 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 1058 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 314 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 1368 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 )

protected

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

Parameters

cov	Matrix which we evaluate Positive Definitiveness

Definition at line 1118 of file ParameterHandlerBase.cpp.

                                                      {
// ********************************************
  if(cov == nullptr){
    cov = &*covMatrix;
    MACH3LOG_WARN("Passed nullptr to cov matrix in {}", matrixName);
  }
 
  M3::MakeMatrixPosDef(cov);
}

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 1520 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 1106 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("============================================================");
}

PrintNominal()

void ParameterHandlerBase::PrintNominal

(

)

const

Print prior value for every parameter.

Definition at line 799 of file ParameterHandlerBase.cpp.

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

PrintNominalCurrProp()

void ParameterHandlerBase::PrintNominalCurrProp

(

)

const

Print prior, current and proposed value for each parameter.

Definition at line 809 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]);
  }
}

PrintParameters()

void ParameterHandlerBase::PrintParameters ( ) const

inline

Warning

only for backward compatibility

Todo:

remove it

Definition at line 253 of file ParameterHandlerBase.h.

{PrintNominalCurrProp();};

ProposeStep()

void ParameterHandlerBase::ProposeStep ( )

virtual

Generate a new proposed state.

Reimplemented in PyParameterHandlerBase.

Definition at line 655 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 593 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 693 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 474 of file ParameterHandlerBase.cpp.

                                                          {
// ********************************************
  _fNames = std::vector<std::string>(SizeVec);
  _fFancyNames = std::vector<std::string>(SizeVec);
  _fPreFitValue = std::vector<double>(SizeVec);
  _fError = std::vector<double>(SizeVec);
  _fCurrVal = std::vector<double>(SizeVec);
  _fPropVal = std::vector<M3::float_t>(SizeVec);
  _fLowBound = std::vector<double>(SizeVec);
  _fUpBound = std::vector<double>(SizeVec);
  _fFlatPrior = std::vector<bool>(SizeVec);
  _fIndivStepScale = std::vector<double>(SizeVec);
  _fSampleNames = std::vector<std::vector<std::string>>(_fNumPar);
 
  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++) {
    _fPreFitValue.at(i) = 1.;
    _fError.at(i) = 1.;
    _fCurrVal.at(i) = 0.;
    _fPropVal.at(i) = 0.;
    _fLowBound.at(i) = -999.99;
    _fUpBound.at(i) = 999.99;
    _fFlatPrior.at(i) = false;
    _fIndivStepScale.at(i) = 1.;
    corr_throw[i] = 0.0;
    randParams[i] = 0.0;
  }
 
  _fGlobalStepScale = 1.0;
}

ResetIndivStepScale()

void ParameterHandlerBase::ResetIndivStepScale

(

)

Adaptive Step Tuning Stuff.

Definition at line 1129 of file ParameterHandlerBase.cpp.

                                               {
// ********************************************
  std::vector<double> stepScales(_fNumPar);
  for (int i = 0; i <_fNumPar; i++) {
    stepScales[i] = 1.;
  }
  _fGlobalStepScale = 1.0;
  SetIndivStepScale(stepScales);
}

RetPointer()

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

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 205 of file ParameterHandlerBase.h.

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

SaveAdaptiveToFile()

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

inline

Save adaptive throw matrix to file.

Definition at line 169 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 1443 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 918 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 452 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 972 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 980 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 987 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 1070 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 994 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 1002 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 1009 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 88 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 1089 of file ParameterHandlerBase.cpp.

                                                                               {
// ********************************************
  if (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 1140 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_INFO("Updating individual step scales for non-adapting parameters to cancel global step scale change.");
  MACH3LOG_INFO("Global step scale initial: {}, current: {}", _fGlobalStepScaleInitial, _fGlobalStepScale);
  PrintIndivStepScale();
}

SetName()

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

inline

Set matrix name.

Definition at line 38 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 180 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 512 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 882 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() != 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 645 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 42 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 56 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 93 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 71 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 636 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 941 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 1193 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

(

TMatrixDSym *

cov

)

Use new throw matrix, used in adaptive MCMC.

Definition at line 1159 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 MakePosDef(throwMatrix);
 
  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 1508 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 670 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 536 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();
}

ToggleFixAllParameters()

void ParameterHandlerBase::ToggleFixAllParameters

(

)

Toggle fixing parameters at prior values.

Definition at line 1016 of file ParameterHandlerBase.cpp.

                                                  {
// ********************************************
  // toggle fix/free all parameters
  if(!pca) for (int i = 0; i < _fNumPar; i++) ToggleFixParameter(i);
  else PCAObj->ToggleFixAllParameters(_fNames);
}

ToggleFixParameter() [1/2]

void ParameterHandlerBase::ToggleFixParameter

(

const int

i

)

Toggle fixing parameter at prior values.

Parameters

i	Parameter index

Definition at line 1024 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

)

Toggle fixing parameter at prior values.

Parameters

name	Name of parameter you want to fix

Definition at line 1046 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 1322 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 1217 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 441 of file ParameterHandlerBase.h.

_fError

std::vector<double> ParameterHandlerBase::_fError

protected

Prior error on the parameter.

Definition at line 445 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 433 of file ParameterHandlerBase.h.

_fFlatPrior

std::vector<bool> ParameterHandlerBase::_fFlatPrior

protected

Whether to apply flat prior or not.

Definition at line 453 of file ParameterHandlerBase.h.

_fGlobalStepScale

double ParameterHandlerBase::_fGlobalStepScale

protected

Global step scale applied to all params in this class.

Definition at line 425 of file ParameterHandlerBase.h.

_fGlobalStepScaleInitial

double ParameterHandlerBase::_fGlobalStepScaleInitial

protected

Backup of _fGlobalStepScale for parameters which are skipped during adaption.

Definition at line 461 of file ParameterHandlerBase.h.

_fIndivStepScale

std::vector<double> ParameterHandlerBase::_fIndivStepScale

protected

Individual step scale used by MCMC algorithm.

Definition at line 451 of file ParameterHandlerBase.h.

_fIndivStepScaleInitial

std::vector<double> ParameterHandlerBase::_fIndivStepScaleInitial

protected

Backup of _fIndivStepScale for parameters which are skipped during adaption.

Definition at line 458 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 447 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 431 of file ParameterHandlerBase.h.

_fNumPar

int ParameterHandlerBase::_fNumPar

protected

Number of systematic parameters.

Definition at line 437 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 439 of file ParameterHandlerBase.h.

_fPropVal

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

protected

Proposed value of the parameter.

Definition at line 443 of file ParameterHandlerBase.h.

_fSampleNames

std::vector<std::vector<std::string> > ParameterHandlerBase::_fSampleNames

protected

Tells to which samples object param should be applied.

Definition at line 455 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 449 of file ParameterHandlerBase.h.

_fYAMLDoc

YAML::Node ParameterHandlerBase::_fYAMLDoc

protected

Stores config describing systematics.

Definition at line 435 of file ParameterHandlerBase.h.

AdaptiveHandler

std::unique_ptr<AdaptiveMCMCHandler> ParameterHandlerBase::AdaptiveHandler

protected

Struct containing information about adaption.

Definition at line 479 of file ParameterHandlerBase.h.

CircularBoundsIndex

std::vector<int> ParameterHandlerBase::CircularBoundsIndex

protected

Indices of parameters with circular bounds.

Definition at line 488 of file ParameterHandlerBase.h.

CircularBoundsValues

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

protected

Circular bounds for each parameter (lower, upper)

Definition at line 490 of file ParameterHandlerBase.h.

corr_throw

double* ParameterHandlerBase::corr_throw

protected

Result of multiplication of Cholesky matrix and randParams.

Definition at line 423 of file ParameterHandlerBase.h.

covMatrix

TMatrixDSym* ParameterHandlerBase::covMatrix

protected

The covariance matrix.

Definition at line 411 of file ParameterHandlerBase.h.

doSpecialStepProposal

bool ParameterHandlerBase::doSpecialStepProposal

protected

Check if any of special step proposal were enabled.

Definition at line 403 of file ParameterHandlerBase.h.

FlipParameterIndex

std::vector<int> ParameterHandlerBase::FlipParameterIndex

protected

Indices of parameters with flip symmetry.

Definition at line 484 of file ParameterHandlerBase.h.

FlipParameterPoint

std::vector<double> ParameterHandlerBase::FlipParameterPoint

protected

Central points around which parameters are flipped.

Definition at line 486 of file ParameterHandlerBase.h.

inputFile

const std::string ParameterHandlerBase::inputFile

protected

The input root file we read in.

Definition at line 406 of file ParameterHandlerBase.h.

invCovMatrix

TMatrixDSym* ParameterHandlerBase::invCovMatrix

protected

The inverse covariance matrix.

Definition at line 413 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 415 of file ParameterHandlerBase.h.

matrixName

std::string ParameterHandlerBase::matrixName

protected

Name of cov matrix.

Definition at line 409 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 464 of file ParameterHandlerBase.h.

pca

bool ParameterHandlerBase::pca

protected

perform PCA or not

Definition at line 472 of file ParameterHandlerBase.h.

PCAObj

std::unique_ptr<PCAHandler> ParameterHandlerBase::PCAObj

protected

Struct containing information about PCA.

Definition at line 477 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 428 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 418 of file ParameterHandlerBase.h.

randParams

double* ParameterHandlerBase::randParams

protected

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

Definition at line 421 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 467 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 469 of file ParameterHandlerBase.h.

Tunes

std::unique_ptr<ParameterTunes> ParameterHandlerBase::Tunes

protected

Struct containing information about adaption.

Definition at line 481 of file ParameterHandlerBase.h.

use_adaptive

bool ParameterHandlerBase::use_adaptive

protected

Are we using AMCMC?

Definition at line 474 of file ParameterHandlerBase.h.

---

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

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