MaCh3  2.2.3
Reference Guide
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
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:
[legend]
Collaboration diagram for ParameterHandlerBase:
[legend]

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 ThrowParProp (const double mag=1.)
 Throw the proposed parameter by mag sigma. Should really just have the user specify this throw by having argument double. More...
 
void ThrowParCurr (const double mag=1.)
 Helper function to throw the current parameter by mag sigma. Can study bias in MCMC with this; put different starting parameters. 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 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 UpdateThrowMatrix (TMatrixDSym *cov)
 
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 double * RetPointer (const int iParam)
 DB Pointer return to param position. More...
 
const std::vector< double > & 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...
 
double 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 [12]. More...
 
void AcceptStep () _noexcept_
 Accepted this step. More...
 
void ToggleFixAllParameters ()
 fix parameters at prior values More...
 
void ToggleFixParameter (const int i)
 fix parameter at prior values More...
 
void ToggleFixParameter (const std::string &name)
 Fix 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...
 
adaptive_mcmc::AdaptiveMCMCHandlerGetAdaptiveHandler () 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...
 
PCAHandlerGetPCAHandler () const
 Get pointer for PCAHandler. More...
 
void MatchMaCh3OutputBranches (TTree *PosteriorFile, std::vector< double > &BranchValues, std::vector< std::string > &BranchNames)
 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 xsec_i, this is currently to make things compatible with the Diagnostic tools. More...
 
std::vector< std::string > _fFancyNames
 Fancy name for example rather than xsec_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< double > _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...
 
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< PCAHandlerPCAObj
 Struct containing information about PCA. More...
 
std::unique_ptr< adaptive_mcmc::AdaptiveMCMCHandlerAdaptiveHandler
 Struct containing information about adaption. More...
 
std::unique_ptr< ParameterTunesTunes
 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.

See also
For more details, visit the Wiki.
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
YAMLFileA vector of strings representing the YAML files used for initialisation of matrix
nameMatrix name
thresholdPCA threshold from 0 to 1. Default is -1 and means no PCA
FirstPCAdparFirst PCA parameter that will be decomposed.
LastPCAdparFirst PCA parameter that will be decomposed.

Definition at line 24 of file ParameterHandlerBase.cpp.

25  : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true) {
26 // ********************************************
27  MACH3LOG_INFO("Constructing instance of ParameterHandler using");
28  doSpecialStepProposal = false;
29  for(unsigned int i = 0; i < YAMLFile.size(); i++)
30  {
31  MACH3LOG_INFO("{}", YAMLFile[i]);
32  }
33  MACH3LOG_INFO("as an input");
34 
35  if (threshold < 0 || threshold >= 1) {
36  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");
37  MACH3LOG_INFO("Please specify a number between 0 and 1");
38  MACH3LOG_INFO("You specified: ");
39  MACH3LOG_INFO("Am instead calling the usual non-PCA constructor...");
40  pca = false;
41  }
42 
43  Init(YAMLFile);
44  // Call the innocent helper function
45  if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
46 }
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:25
bool doSpecialStepProposal
Check if any of special step proposal were enabled.
void Init(const std::string &name, const std::string &file)
Initialisation of the class using matrix from root file.
std::string matrixName
Name of cov matrix.
void ConstructPCA(const double eigen_threshold, int FirstPCAdpar, int LastPCAdpar)
CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold.
const std::string inputFile
The input root file we read in.
bool pca
perform PCA or not

◆ 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
nameMatrix name
filePath to matrix root file

Definition at line 5 of file ParameterHandlerBase.cpp.

6  : inputFile(file), pca(false) {
7 // ********************************************
8  MACH3LOG_DEBUG("Constructing instance of ParameterHandler");
9  doSpecialStepProposal = false;
10  if (threshold < 0 || threshold >= 1) {
11  MACH3LOG_INFO("NOTE: {} {}", name, file);
12  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");
13  MACH3LOG_INFO("Please specify a number between 0 and 1");
14  MACH3LOG_INFO("You specified: ");
15  MACH3LOG_INFO("Am instead calling the usual non-PCA constructor...");
16  pca = false;
17  }
18  Init(name, file);
19 
20  // Call the innocent helper function
21  if (pca) ConstructPCA(threshold, FirstPCA, LastPCA);
22 }
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:24

◆ ~ParameterHandlerBase()

ParameterHandlerBase::~ParameterHandlerBase ( )
virtual

Destructor.

Definition at line 50 of file ParameterHandlerBase.cpp.

50  {
51 // ********************************************
52  delete[] randParams;
53  delete[] corr_throw;
54 
55  if (covMatrix != nullptr) delete covMatrix;
56  if (invCovMatrix != nullptr) delete invCovMatrix;
57  if (throwMatrix != nullptr) delete throwMatrix;
58  for(int i = 0; i < _fNumPar; i++) {
59  delete[] throwMatrixCholDecomp[i];
60  }
61  delete[] throwMatrixCholDecomp;
62 }
double * randParams
Random number taken from gaussian around prior error used for corr_throw.
double ** throwMatrixCholDecomp
Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss.
TMatrixDSym * invCovMatrix
The inverse covariance matrix.
TMatrixDSym * throwMatrix
Matrix which we use for step proposal before Cholesky decomposition (not actually used for step propo...
int _fNumPar
Number of systematic parameters.
TMatrixDSym * covMatrix
The covariance matrix.
double * corr_throw
Result of multiplication of Cholesky matrix and randParams.

Member Function Documentation

◆ AcceptStep()

void ParameterHandlerBase::AcceptStep ( )

Accepted this step.

Definition at line 656 of file ParameterHandlerBase.cpp.

656  {
657 // ********************************************
658  if (!pca) {
659  #ifdef MULTITHREAD
660  #pragma omp parallel for
661  #endif
662  for (int i = 0; i < _fNumPar; ++i) {
663  // Update state so that current state is proposed state
664  _fCurrVal[i] = _fPropVal[i];
665  }
666  } else {
667  PCAObj->AcceptStep();
668  }
669 }
std::unique_ptr< PCAHandler > PCAObj
Struct containing information about PCA.
std::vector< double > _fCurrVal
Current value of the parameter.
std::vector< double > _fPropVal
Proposed value of the parameter.

◆ AppliesToSample()

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

Check if parameter is affecting given sample name.

Parameters
SystIndexnumber of parameter
SampleNameThe Sample name used to filter parameters.

Definition at line 1235 of file ParameterHandlerBase.cpp.

1235  {
1236 // ********************************************
1237  // Empty means apply to all
1238  if (_fSampleNames[SystIndex].size() == 0) return true;
1239 
1240  // Make a copy and to lower case to not be case sensitive
1241  std::string SampleNameCopy = SampleName;
1242  std::transform(SampleNameCopy.begin(), SampleNameCopy.end(), SampleNameCopy.begin(), ::tolower);
1243 
1244  // Check for unsupported wildcards in SampleNameCopy
1245  if (SampleNameCopy.find('*') != std::string::npos) {
1246  MACH3LOG_ERROR("Wildcards ('*') are not supported in sample name: '{}'", SampleName);
1247  throw MaCh3Exception(__FILE__ , __LINE__ );
1248  }
1249 
1250  bool Applies = false;
1251 
1252  for (size_t i = 0; i < _fSampleNames[SystIndex].size(); i++) {
1253  // Convert to low case to not be case sensitive
1254  std::string pattern = _fSampleNames[SystIndex][i];
1255  std::transform(pattern.begin(), pattern.end(), pattern.begin(), ::tolower);
1256 
1257  // Replace '*' in the pattern with '.*' for regex matching
1258  std::string regexPattern = "^" + std::regex_replace(pattern, std::regex("\\*"), ".*") + "$";
1259  try {
1260  std::regex regex(regexPattern);
1261  if (std::regex_match(SampleNameCopy, regex)) {
1262  Applies = true;
1263  break;
1264  }
1265  } catch (const std::regex_error& e) {
1266  // Handle regex error (for invalid patterns)
1267  MACH3LOG_ERROR("Regex error: {}", e.what());
1268  }
1269  }
1270  return Applies;
1271 }
int size
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:27
Custom exception class for MaCh3 errors.
std::vector< std::vector< std::string > > _fSampleNames
Tells to which samples object param should be applied.

◆ CalcLikelihood()

double ParameterHandlerBase::CalcLikelihood ( ) const

Calc penalty term based on inverted covariance matrix.

Definition at line 758 of file ParameterHandlerBase.cpp.

758  {
759 // ********************************************
760  double logL = 0.0;
761  #ifdef MULTITHREAD
762  #pragma omp parallel for reduction(+:logL)
763  #endif
764  for(int i = 0; i < _fNumPar; ++i){
765  if(_fFlatPrior[i]){
766  //HW: Flat prior, no need to calculate anything
767  continue;
768  }
769 
770  #ifdef MULTITHREAD
771  #pragma omp simd
772  #endif
773  for (int j = 0; j <= i; ++j) {
774  if (!_fFlatPrior[j]) {
775  //KS: Since matrix is symmetric we can calculate non diagonal elements only once and multiply by 2, can bring up to factor speed decrease.
776  double scale = (i != j) ? 1. : 0.5;
777  logL += scale * (_fPropVal[i] - _fPreFitValue[i])*(_fPropVal[j] - _fPreFitValue[j])*InvertCovMatrix[i][j];
778  }
779  }
780  }
781  return logL;
782 }
std::vector< bool > _fFlatPrior
Whether to apply flat prior or not.
std::vector< std::vector< double > > InvertCovMatrix
KS: Same as above but much faster as TMatrixDSym cache miss.
std::vector< double > _fPreFitValue
Parameter value dictated by the prior model. Based on it penalty term is calculated.

◆ CheckBounds()

int ParameterHandlerBase::CheckBounds ( ) const

Check if parameters were proposed outside physical boundary.

Definition at line 785 of file ParameterHandlerBase.cpp.

785  {
786 // ********************************************
787  int NOutside = 0;
788  #ifdef MULTITHREAD
789  #pragma omp parallel for reduction(+:NOutside)
790  #endif
791  for (int i = 0; i < _fNumPar; ++i){
792  if(_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]){
793  NOutside++;
794  }
795  }
796  return NOutside;
797 }
std::vector< double > _fLowBound
Lowest physical bound, parameter will not be able to go beyond it.
std::vector< double > _fUpBound
Upper physical bound, parameter will not be able to go beyond it.

◆ 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 673 of file ParameterHandlerBase.cpp.

673  {
674 // *************************************
675  if(_fPropVal[index] > UpBound) {
676  _fPropVal[index] = LowBound + std::fmod(_fPropVal[index] - UpBound, UpBound - LowBound);
677  } else if (_fPropVal[index] < LowBound) {
678  _fPropVal[index] = UpBound - std::fmod(LowBound - _fPropVal[index], UpBound - LowBound);
679  }
680 }

◆ 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_thresholdPCA threshold from 0 to 1. Default is -1 and means no PCA
FirstPCAdparFirst PCA parameter that will be decomposed.
LastPCAdparFirst PCA parameter that will be decomposed.
See also
For more details, visit the Wiki.

Definition at line 65 of file ParameterHandlerBase.cpp.

65  {
66 // ********************************************
67  if(AdaptiveHandler) {
68  MACH3LOG_ERROR("Adaption has been enabled and now trying to enable PCA. Right now both configuration don't work with each other");
69  throw MaCh3Exception(__FILE__ , __LINE__ );
70  }
71 
72  PCAObj = std::make_unique<PCAHandler>();
73  //Check whether first and last pcadpar are set and if not just PCA everything
74  if(FirstPCAdpar == -999 || LastPCAdpar == -999){
75  if(FirstPCAdpar == -999 && LastPCAdpar == -999){
76  FirstPCAdpar = 0;
77  LastPCAdpar = covMatrix->GetNrows()-1;
78  }
79  else{
80  MACH3LOG_ERROR("You must either leave FirstPCAdpar and LastPCAdpar at -999 or set them both to something");
81  throw MaCh3Exception(__FILE__ , __LINE__ );
82  }
83  }
84 
85  PCAObj->ConstructPCA(covMatrix, FirstPCAdpar, LastPCAdpar, eigen_threshold, _fNumPar);
86  PCAObj->SetupPointers(&_fCurrVal, &_fPropVal);
87  // Make a note that we have now done PCA
88  pca = true;
89 }
std::unique_ptr< adaptive_mcmc::AdaptiveMCMCHandler > AdaptiveHandler
Struct containing information about adaption.

◆ CorrelateSteps()

void ParameterHandlerBase::CorrelateSteps ( )

Use Cholesky throw matrix for better step proposal.

Definition at line 634 of file ParameterHandlerBase.cpp.

634  {
635 // ************************************************
636  //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
638 
639  // If not doing PCA
640  if (!pca) {
641  #ifdef MULTITHREAD
642  #pragma omp parallel for
643  #endif
644  for (int i = 0; i < _fNumPar; ++i) {
645  if (!IsParameterFixed(i) > 0.) {
647  }
648  }
649  // If doing PCA throw uncorrelated in PCA basis (orthogonal basis by definition)
650  } else {
652  }
653 }
double _fGlobalStepScale
Global step scale applied to all params in this class.
bool IsParameterFixed(const int i) const
Is parameter fixed or not.
std::vector< double > _fIndivStepScale
Individual step scale used by MCMC algorithm.
void MatrixVectorMulti(double *_restrict_ VecMulti, double **_restrict_ matrix, const double *_restrict_ vector, const int n)
KS: Custom function to perform multiplication of matrix and vector with multithreading.

◆ EnableSpecialProposal()

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

Enable special proposal.

Definition at line 296 of file ParameterHandlerBase.cpp.

296  {
297 // ********************************************
298  doSpecialStepProposal = true;
299 
300  bool CircEnabled = false;
301  bool FlipEnabled = false;
302 
303  if (param["CircularBounds"]) {
304  CircEnabled = true;
305  }
306 
307  if (param["FlipParameter"]) {
308  FlipEnabled = true;
309  }
310 
311  if (!CircEnabled && !FlipEnabled) {
312  MACH3LOG_ERROR("None of Special Proposal were enabled even though param {}, has SpecialProposal entry in Yaml", GetParFancyName(Index));
313  throw MaCh3Exception(__FILE__, __LINE__);
314  }
315 
316  if (CircEnabled) {
317  CircularBoundsIndex.push_back(Index);
318  CircularBoundsValues.push_back(Get<std::pair<double, double>>(param["CircularBounds"], __FILE__, __LINE__));
319  MACH3LOG_INFO("Enabling CircularBounds for parameter {} with range [{}, {}]",
320  GetParFancyName(Index),
321  CircularBoundsValues.back().first,
322  CircularBoundsValues.back().second);
323  // KS: Make sure circular bounds are within physical bounds. If we are outside of physics bound MCMC will never explore such phase space region
324  if (CircularBoundsValues.back().first < _fLowBound.at(Index) || CircularBoundsValues.back().second > _fUpBound.at(Index)) {
325  MACH3LOG_ERROR("Circular bounds [{}, {}] for parameter {} exceed physical bounds [{}, {}]",
326  CircularBoundsValues.back().first, CircularBoundsValues.back().second,
327  GetParFancyName(Index),
328  _fLowBound.at(Index), _fUpBound.at(Index));
329  throw MaCh3Exception(__FILE__, __LINE__);
330  }
331  }
332 
333  if (FlipEnabled) {
334  FlipParameterIndex.push_back(Index);
335  FlipParameterPoint.push_back(Get<double>(param["FlipParameter"], __FILE__, __LINE__));
336  MACH3LOG_INFO("Enabling Flipping for parameter {} with value {}",
337  GetParFancyName(Index),
338  FlipParameterPoint.back());
339  }
340 
341  if (CircEnabled && FlipEnabled) {
342  if (FlipParameterPoint.back() < CircularBoundsValues.back().first || FlipParameterPoint.back() > CircularBoundsValues.back().second) {
343  MACH3LOG_ERROR("FlipParameter value {} for parameter {} is outside the CircularBounds [{}, {}]",
344  FlipParameterPoint.back(), GetParFancyName(Index), CircularBoundsValues.back().first, CircularBoundsValues.back().second);
345  throw MaCh3Exception(__FILE__, __LINE__);
346  }
347 
348  const double low = CircularBoundsValues.back().first;
349  const double high = CircularBoundsValues.back().second;
350 
351  // Sanity check: ensure flipping any x in [low, high] keeps the result in [low, high]
352  const double flipped_low = 2 * FlipParameterPoint.back() - low;
353  const double flipped_high = 2 * FlipParameterPoint.back() - high;
354  const double min_flip = std::min(flipped_low, flipped_high);
355  const double max_flip = std::max(flipped_low, flipped_high);
356 
357  if (min_flip < low || max_flip > high) {
358  MACH3LOG_ERROR("Flipping about point {} for parameter {} would leave circular bounds [{}, {}]",
359  FlipParameterPoint.back(), GetParFancyName(Index), low, high);
360  throw MaCh3Exception(__FILE__, __LINE__);
361  }
362  }
363 }
Type Get(const YAML::Node &node, const std::string File, const int Line)
Get content of config file.
Definition: YamlHelper.h:266
std::vector< int > CircularBoundsIndex
Indices of parameters with circular bounds.
std::vector< int > FlipParameterIndex
Indices of parameters with flip symmetry.
std::vector< double > FlipParameterPoint
Central points around which parameters are flipped.
std::vector< std::pair< double, double > > CircularBoundsValues
Circular bounds for each parameter (lower, upper)
std::string GetParFancyName(const int i) const
Get fancy name of the Parameter.

◆ FlipParameterValue()

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

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

Parameters
indexparameter index you want to flip
FlipPointValue around which flipping is done

Definition at line 683 of file ParameterHandlerBase.cpp.

683  {
684 // *************************************
685  if(random_number[0]->Uniform() < 0.5) {
686  _fPropVal[index] = 2 * FlipPoint - _fPropVal[index];
687  }
688 }
std::vector< std::unique_ptr< TRandom3 > > random_number
KS: Set Random numbers for each thread so each thread has different seed.

◆ GetParPropVec()

const std::vector<double>& 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 237 of file ParameterHandlerBase.h.

237 {return _fPropVal;}

◆ GetPCAHandler()

PCAHandler* ParameterHandlerBase::GetPCAHandler ( ) const
inline

Get pointer for PCAHandler.

Definition at line 358 of file ParameterHandlerBase.h.

358  {
359  if (!pca) {
360  MACH3LOG_ERROR("Am not running in PCA mode");
361  throw MaCh3Exception(__FILE__ , __LINE__ );
362  }
363  return PCAObj.get();
364  }

◆ 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 92 of file ParameterHandlerBase.cpp.

92  {
93 // ********************************************
94  // Set the covariance matrix from input ROOT file (e.g. flux, ND280, NIWG)
95  TFile *infile = new TFile(file.c_str(), "READ");
96  if (infile->IsZombie()) {
97  MACH3LOG_ERROR("Could not open input covariance ROOT file {} !!!", file);
98  MACH3LOG_ERROR("Was about to retrieve matrix with name {}", name);
99  throw MaCh3Exception(__FILE__ , __LINE__ );
100  }
101 
102  TMatrixDSym *CovMat = static_cast<TMatrixDSym*>(infile->Get(name.c_str()));
103 
104  if (!CovMat) {
105  MACH3LOG_ERROR("Could not find covariance matrix name {} in file {}", name, file);
106  MACH3LOG_ERROR("Are you really sure {} exists in the file?", name);
107  throw MaCh3Exception(__FILE__ , __LINE__ );
108  }
109 
110  PrintLength = 35;
111 
112  const int nThreads = M3::GetNThreads();
113  //KS: set Random numbers for each thread so each thread has different seed
114  //or for one thread if without MULTITHREAD
115  random_number.reserve(nThreads);
116  for (int iThread = 0; iThread < nThreads; iThread++) {
117  random_number.emplace_back(std::make_unique<TRandom3>(0));
118  }
119  // Not using adaptive by default
120  use_adaptive = false;
121  // Set the covariance matrix
122  _fNumPar = CovMat->GetNrows();
123 
124  InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
125  throwMatrixCholDecomp = new double*[_fNumPar]();
126  // Set the defaults to true
127  for(int i = 0; i < _fNumPar; i++) {
128  throwMatrixCholDecomp[i] = new double[_fNumPar]();
129  for (int j = 0; j < _fNumPar; j++) {
130  throwMatrixCholDecomp[i][j] = 0.;
131  }
132  }
133  SetName(name);
134  MakePosDef(CovMat);
135  SetCovMatrix(CovMat);
136  if (_fNumPar <= 0) {
137  MACH3LOG_CRITICAL("Covariance matrix {} has {} entries!", GetName(), _fNumPar);
138  throw MaCh3Exception(__FILE__ , __LINE__ );
139  }
140 
142 
143  infile->Close();
144 
145  MACH3LOG_INFO("Created covariance matrix named: {}", GetName());
146  MACH3LOG_INFO("from file: {}", file);
147  delete infile;
148 }
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:28
bool use_adaptive
Are we using AMCMC?
void ReserveMemory(const int size)
Initialise vectors with parameters information.
void MakePosDef(TMatrixDSym *cov=nullptr)
Make matrix positive definite by adding small values to diagonal, necessary for inverting matrix.
int PrintLength
KS: This is used when printing parameters, sometimes we have super long parameters name,...
std::string GetName() const
Get name of covariance.
void SetName(const std::string &name)
Set matrix name.
void SetCovMatrix(TMatrixDSym *cov)
Set covariance matrix.
int GetNThreads()
number of threads which we need for example for TRandom3
Definition: Monitor.cpp:357

◆ Init() [2/2]

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

Initialisation of the class using config.

Parameters
YAMLFileA vector of strings representing the YAML files used for initialisation of matrix

Definition at line 153 of file ParameterHandlerBase.cpp.

153  {
154 // ********************************************
155  _fYAMLDoc["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
156  for(unsigned int i = 0; i < YAMLFile.size(); i++)
157  {
158  YAML::Node YAMLDocTemp = M3OpenConfig(YAMLFile[i]);
159  for (const auto& item : YAMLDocTemp["Systematics"]) {
160  _fYAMLDoc["Systematics"].push_back(item);
161  }
162  }
163 
164  const int nThreads = M3::GetNThreads();
165  //KS: set Random numbers for each thread so each thread has different seed
166  //or for one thread if without MULTITHREAD
167  random_number.reserve(nThreads);
168  for (int iThread = 0; iThread < nThreads; iThread++) {
169  random_number.emplace_back(std::make_unique<TRandom3>(0));
170  }
171  PrintLength = 35;
172 
173  // Set the covariance matrix
174  _fNumPar = int(_fYAMLDoc["Systematics"].size());
175 
176  use_adaptive = false;
177 
178  InvertCovMatrix.resize(_fNumPar, std::vector<double>(_fNumPar, 0.0));
179  throwMatrixCholDecomp = new double*[_fNumPar]();
180  for(int i = 0; i < _fNumPar; i++) {
181  throwMatrixCholDecomp[i] = new double[_fNumPar]();
182  for (int j = 0; j < _fNumPar; j++) {
183  throwMatrixCholDecomp[i][j] = 0.;
184  }
185  }
187 
188  TMatrixDSym* _fCovMatrix = new TMatrixDSym(_fNumPar);
189  int i = 0;
190  std::vector<std::map<std::string,double>> Correlations(_fNumPar);
191  std::map<std::string, int> CorrNamesMap;
192 
193  //ETA - read in the systematics. Would be good to add in some checks to make sure
194  //that there are the correct number of entries i.e. are the _fNumPar for Names,
195  //PreFitValues etc etc.
196  for (auto const &param : _fYAMLDoc["Systematics"])
197  {
198  _fFancyNames[i] = Get<std::string>(param["Systematic"]["Names"]["FancyName"], __FILE__ , __LINE__);
199  _fPreFitValue[i] = Get<double>(param["Systematic"]["ParameterValues"]["PreFitValue"], __FILE__ , __LINE__);
200  _fIndivStepScale[i] = Get<double>(param["Systematic"]["StepScale"]["MCMC"], __FILE__ , __LINE__);
201  _fError[i] = Get<double>(param["Systematic"]["Error"], __FILE__ , __LINE__);
202  _fSampleNames[i] = GetFromManager<std::vector<std::string>>(param["Systematic"]["SampleNames"], {}, __FILE__, __LINE__);
203  if(_fError[i] <= 0) {
204  MACH3LOG_ERROR("Error for param {}({}) is negative and equal to {}", _fFancyNames[i], i, _fError[i]);
205  throw MaCh3Exception(__FILE__ , __LINE__ );
206  }
207  //ETA - a bit of a fudge but works
208  auto TempBoundsVec = GetBounds(param["Systematic"]["ParameterBounds"]);
209  _fLowBound[i] = TempBoundsVec[0];
210  _fUpBound[i] = TempBoundsVec[1];
211 
212  //ETA - now for parameters which are optional and have default values
213  _fFlatPrior[i] = GetFromManager<bool>(param["Systematic"]["FlatPrior"], false, __FILE__ , __LINE__);
214 
215  // 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
216  if(GetFromManager<bool>(param["Systematic"]["FixParam"], false, __FILE__ , __LINE__)) {
218  }
219 
220  if(param["Systematic"]["SpecialProposal"]) {
221  EnableSpecialProposal(param["Systematic"]["SpecialProposal"], i);
222  }
223 
224  //Fill the map to get the correlations later as well
225  CorrNamesMap[param["Systematic"]["Names"]["FancyName"].as<std::string>()]=i;
226 
227  //Also loop through the correlations
228  if(param["Systematic"]["Correlations"]) {
229  for(unsigned int Corr_i = 0; Corr_i < param["Systematic"]["Correlations"].size(); ++Corr_i){
230  for (YAML::const_iterator it = param["Systematic"]["Correlations"][Corr_i].begin(); it!=param["Systematic"]["Correlations"][Corr_i].end();++it) {
231  Correlations[i][it->first.as<std::string>()] = it->second.as<double>();
232  }
233  }
234  }
235  i++;
236  } // end loop over para
237  if(i != _fNumPar) {
238  MACH3LOG_CRITICAL("Inconsistent number of params in Yaml {} vs {}, this indicate wrong syntax", i, i, _fNumPar);
239  throw MaCh3Exception(__FILE__ , __LINE__ );
240  }
241  // ETA Now that we've been through all systematic let's fill the covmatrix
242  //This makes the root TCov from YAML
243  for(int j = 0; j < _fNumPar; j++) {
244  (*_fCovMatrix)(j, j) = _fError[j]*_fError[j];
245  //Get the map of parameter name to correlation from the Correlations object
246  for (auto const& pair : Correlations[j]) {
247  auto const& key = pair.first;
248  auto const& val = pair.second;
249  int index = -1;
250  //If you found the parameter name then get the index
251  if (CorrNamesMap.find(key) != CorrNamesMap.end()) {
252  index = CorrNamesMap[key];
253  } else {
254  MACH3LOG_ERROR("Parameter {} not in list! Check your spelling?", key);
255  throw MaCh3Exception(__FILE__ , __LINE__ );
256  }
257  double Corr1 = val;
258  double Corr2 = 0;
259  if(Correlations[index].find(_fFancyNames[j]) != Correlations[index].end()) {
260  Corr2 = Correlations[index][_fFancyNames[j]];
261  //Do they agree to better than float precision?
262  if(std::abs(Corr2 - Corr1) > FLT_EPSILON) {
263  MACH3LOG_ERROR("Correlations are not equal between {} and {}", _fFancyNames[j], key);
264  MACH3LOG_ERROR("Got : {} and {}", Corr2, Corr1);
265  throw MaCh3Exception(__FILE__ , __LINE__ );
266  }
267  } else {
268  MACH3LOG_ERROR("Correlation does not appear reciprocally between {} and {}", _fFancyNames[j], key);
269  throw MaCh3Exception(__FILE__ , __LINE__ );
270  }
271  (*_fCovMatrix)(j, index)= (*_fCovMatrix)(index, j) = Corr1*_fError[j]*_fError[index];
272  }
273  }
274 
275  //Now make positive definite
276  MakePosDef(_fCovMatrix);
277  SetCovMatrix(_fCovMatrix);
278 
279  if (_fNumPar <= 0) {
280  MACH3LOG_ERROR("ParameterHandler object has {} systematics!", _fNumPar);
281  throw MaCh3Exception(__FILE__ , __LINE__ );
282  }
283 
284  Tunes = std::make_unique<ParameterTunes>(_fYAMLDoc["Systematics"]);
285 
286  MACH3LOG_INFO("Created covariance matrix from files: ");
287  for(const auto &file : YAMLFile){
288  MACH3LOG_INFO("{} ", file);
289  }
290  MACH3LOG_INFO("----------------");
291  MACH3LOG_INFO("Found {} systematics parameters in total", _fNumPar);
292  MACH3LOG_INFO("----------------");
293 }
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:561
#define GetBounds(filename)
Definition: YamlHelper.h:562
void ToggleFixParameter(const int i)
fix parameter at prior values
std::unique_ptr< ParameterTunes > Tunes
Struct containing information about adaption.
std::vector< std::string > _fFancyNames
Fancy name for example rather than xsec_0 it is MAQE, useful for human reading.
std::vector< double > _fError
Prior error on the parameter.
YAML::Node _fYAMLDoc
Stores config describing systematics.
void EnableSpecialProposal(const YAML::Node &param, const int Index)
Enable special proposal.

◆ InitialiseAdaption()

void ParameterHandlerBase::InitialiseAdaption ( const YAML::Node &  adapt_manager)

Initialise adaptive MCMC.

Parameters
adapt_managerNode having from which we load all adaptation options

Definition at line 1062 of file ParameterHandlerBase.cpp.

1062  {
1063 // ********************************************
1064  if(PCAObj){
1065  MACH3LOG_ERROR("PCA has been enabled and now trying to enable Adaption. Right now both configuration don't work with each other");
1066  throw MaCh3Exception(__FILE__ , __LINE__ );
1067  }
1068  if(AdaptiveHandler){
1069  MACH3LOG_ERROR("Adaptive Handler has already been initialise can't do it again so skipping.");
1070  return;
1071  }
1072  AdaptiveHandler = std::make_unique<adaptive_mcmc::AdaptiveMCMCHandler>();
1073  // Now we read the general settings [these SHOULD be common across all matrices!]
1074  bool success = AdaptiveHandler->InitFromConfig(adapt_manager, matrixName, &_fCurrVal, &_fError);
1075  if(!success) return;
1076  AdaptiveHandler->Print();
1077 
1078  // Next let"s check for external matrices
1079  // We"re going to grab this info from the YAML manager
1080  if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["UseExternalMatrix"], false, __FILE__ , __LINE__)) {
1081  MACH3LOG_WARN("Not using external matrix for {}, initialising adaption from scratch", matrixName);
1082  // If we don't have a covariance matrix to start from for adaptive tune we need to make one!
1083  use_adaptive = true;
1084  AdaptiveHandler->CreateNewAdaptiveCovariance();
1085  return;
1086  }
1087 
1088  // Finally, we accept that we want to read the matrix from a file!
1089  auto external_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixFileName"], "", __FILE__ , __LINE__);
1090  auto external_matrix_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixName"], "", __FILE__ , __LINE__);
1091  auto external_mean_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMeansName"], "", __FILE__ , __LINE__);
1092 
1093  AdaptiveHandler->SetThrowMatrixFromFile(external_file_name, external_matrix_name, external_mean_name, use_adaptive);
1094  SetThrowMatrix(AdaptiveHandler->GetAdaptiveCovariance());
1096 
1097  MACH3LOG_INFO("Successfully Set External Throw Matrix Stored in {}", external_file_name);
1098 }
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:26
void ResetIndivStepScale()
Adaptive Step Tuning Stuff.
void SetThrowMatrix(TMatrixDSym *cov)
Use new throw matrix, used in adaptive MCMC.

◆ IsParameterFixed() [1/2]

bool ParameterHandlerBase::IsParameterFixed ( const int  i) const
inline

Is parameter fixed or not.

Parameters
iParameter index

Definition at line 321 of file ParameterHandlerBase.h.

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

◆ IsParameterFixed() [2/2]

bool ParameterHandlerBase::IsParameterFixed ( const std::string &  name) const

Is parameter fixed or not.

Parameters
nameName of parameter you want to check if is fixed

Definition at line 936 of file ParameterHandlerBase.cpp.

936  {
937 // ********************************************
938  const int Index = GetParIndex(name);
939  if(Index != M3::_BAD_INT_) {
940  return IsParameterFixed(Index);
941  }
942 
943  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore don't know if it fixed", name);
944  return false;
945 }
int GetParIndex(const std::string &name) const
Get index based on name.
constexpr static const int _BAD_INT_
Default value used for int initialisation.
Definition: Core.h:48

◆ IsPCA()

bool ParameterHandlerBase::IsPCA ( ) const
inline

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

Definition at line 337 of file ParameterHandlerBase.h.

337 { 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 1135 of file ParameterHandlerBase.cpp.

1135  {
1136 // ********************************************
1137  // Want to get cov' = (cov_sym+cov_polar)/2
1138  // cov_sym=(cov+cov^T)/2
1139  // cov_polar-> SVD cov to cov=USV^T then cov_polar=VSV^T
1140 
1141  //Get frob norm of cov
1142  // Double_t cov_norm=cov->E2Norm();
1143 
1144  TMatrixDSym* cov_trans = cov;
1145  cov_trans->T();
1146  TMatrixDSym cov_sym = 0.5*(*cov+*cov_trans); //If cov is symmetric does nothing, otherwise just ensures symmetry
1147 
1148  //Do SVD to get polar form
1149  TDecompSVD cov_sym_svd=TDecompSVD(cov_sym);
1150  if(!cov_sym_svd.Decompose()){
1151  MACH3LOG_ERROR("Cannot do SVD on input matrix!");
1152  throw MaCh3Exception(__FILE__ , __LINE__ );
1153  }
1154 
1155  TMatrixD cov_sym_v = cov_sym_svd.GetV();
1156  TMatrixD cov_sym_vt = cov_sym_v;
1157  cov_sym_vt.T();
1158  //SVD returns as vector (grrr) so need to get into matrix form for multiplying!
1159  TVectorD cov_sym_sigvect = cov_sym_svd.GetSig();
1160 
1161  const Int_t nCols = cov_sym_v.GetNcols(); //square so only need rows hence lack of cols
1162  TMatrixDSym cov_sym_sig(nCols);
1163  TMatrixDDiag cov_sym_sig_diag(cov_sym_sig);
1164  cov_sym_sig_diag=cov_sym_sigvect;
1165 
1166  //Can finally get H=VSV
1167  TMatrixDSym cov_sym_polar = cov_sym_sig.SimilarityT(cov_sym_vt);//V*S*V^T (this took forver to find!)
1168 
1169  //Now we can construct closest approximater Ahat=0.5*(B+H)
1170  TMatrixDSym cov_closest_approx = 0.5*(cov_sym+cov_sym_polar);//Not fully sure why this is even needed since symmetric B -> U=V
1171  //Get norm of transformed
1172  // Double_t approx_norm=cov_closest_approx.E2Norm();
1173  //MACH3LOG_INFO("Initial Norm: {:.6f} | Norm after transformation: {:.6f} | Ratio: {:.6f}", cov_norm, approx_norm, cov_norm / approx_norm);
1174 
1175  *cov = cov_closest_approx;
1176  //Now can just add a makeposdef!
1177  MakePosDef(cov);
1178 }

◆ MakePosDef()

void ParameterHandlerBase::MakePosDef ( TMatrixDSym *  cov = nullptr)
protected

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

Parameters
covMatrix which we evaluate Positive Definitiveness

Definition at line 996 of file ParameterHandlerBase.cpp.

996  {
997 // ********************************************
998  if(cov == nullptr){
999  cov = &*covMatrix;
1000  MACH3LOG_WARN("Passed nullptr to cov matrix in {}", matrixName);
1001  }
1002 
1003  M3::MakeMatrixPosDef(cov);
1004 }
void MakeMatrixPosDef(TMatrixDSym *cov)
Makes sure that matrix is positive-definite by adding a small number to on-diagonal elements.

◆ MatchMaCh3OutputBranches()

void ParameterHandlerBase::MatchMaCh3OutputBranches ( TTree *  PosteriorFile,
std::vector< double > &  BranchValues,
std::vector< std::string > &  BranchNames 
)

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

Parameters
PosteriorFilePointer to the ROOT TTree from MaCh3 fit.
[out]BranchValuesVector to store the values of the branches (resized inside).
[out]BranchNamesVector to store the names of the branches (resized inside).
Exceptions
MaCh3Exceptionif any parameter branch is uninitialized.
Parameters
PosteriorFilePointer to the ROOT TTree from MaCh3 fit.
SystematicPointer to the systematic parameter handler.
[out]BranchValuesVector to store the values of the branches (resized inside).
[out]BranchNamesVector to store the names of the branches (resized inside).
Exceptions
MaCh3Exceptionif any parameter branch is uninitialized.

Definition at line 1295 of file ParameterHandlerBase.cpp.

1297  {
1298 // *************************************
1299  BranchValues.resize(GetNumParams());
1300  BranchNames.resize(GetNumParams());
1301 
1302  for (int i = 0; i < GetNumParams(); ++i) {
1303  BranchNames[i] = GetParName(i);
1304  if (!PosteriorFile->GetBranch(BranchNames[i].c_str())) {
1305  MACH3LOG_ERROR("Branch '{}' does not exist in the TTree!", BranchNames[i]);
1306  throw MaCh3Exception(__FILE__, __LINE__);
1307  }
1308  PosteriorFile->SetBranchStatus(BranchNames[i].c_str(), true);
1309  PosteriorFile->SetBranchAddress(BranchNames[i].c_str(), &BranchValues[i]);
1310  }
1311 }
std::vector< TString > BranchNames
int GetNumParams() const
Get total number of parameters.
std::string GetParName(const int i) const
Get name of parameter.

◆ PrintIndivStepScale()

void ParameterHandlerBase::PrintIndivStepScale ( ) const

Print step scale for each parameter.

Definition at line 984 of file ParameterHandlerBase.cpp.

984  {
985 // ********************************************
986  MACH3LOG_INFO("============================================================");
987  MACH3LOG_INFO("{:<{}} | {:<11}", "Parameter:", PrintLength, "Step scale:");
988  for (int iParam = 0; iParam < _fNumPar; iParam++) {
989  MACH3LOG_INFO("{:<{}} | {:<11}", _fFancyNames[iParam].c_str(), PrintLength, _fIndivStepScale[iParam]);
990  }
991  MACH3LOG_INFO("============================================================");
992 }

◆ PrintNominal()

void ParameterHandlerBase::PrintNominal ( ) const

Print prior value for every parameter.

Definition at line 730 of file ParameterHandlerBase.cpp.

730  {
731 // ********************************************
732  MACH3LOG_INFO("Prior values for {} ParameterHandler:", GetName());
733  for (int i = 0; i < _fNumPar; i++) {
734  MACH3LOG_INFO(" {} {} ", GetParFancyName(i), GetParInit(i));
735  }
736 }
double GetParInit(const int i) const
Get prior parameter value.

◆ PrintNominalCurrProp()

void ParameterHandlerBase::PrintNominalCurrProp ( ) const

Print prior, current and proposed value for each parameter.

Definition at line 740 of file ParameterHandlerBase.cpp.

740  {
741 // ********************************************
742  MACH3LOG_INFO("Printing parameters for {}", GetName());
743  // Dump out the PCA parameters too
744  if (pca) {
745  PCAObj->Print();
746  }
747  MACH3LOG_INFO("{:<30} {:<10} {:<10} {:<10}", "Name", "Prior", "Current", "Proposed");
748  for (int i = 0; i < _fNumPar; ++i) {
749  MACH3LOG_INFO("{:<30} {:<10.2f} {:<10.2f} {:<10.2f}", GetParFancyName(i), _fPreFitValue[i], _fCurrVal[i], _fPropVal[i]);
750  }
751 }

◆ PrintParameters()

void ParameterHandlerBase::PrintParameters ( ) const
inline
Warning
only for backward compatibility
Todo:
remove it

Definition at line 294 of file ParameterHandlerBase.h.

void PrintNominalCurrProp() const
Print prior, current and proposed value for each parameter.

◆ ProposeStep()

void ParameterHandlerBase::ProposeStep ( )
virtual

Generate a new proposed state.

Reimplemented in PyParameterHandlerBase.

Definition at line 560 of file ParameterHandlerBase.cpp.

560  {
561 // ************************************************
562  // Make the random numbers for the step proposal
563  Randomize();
564  CorrelateSteps();
565 
566  // KS: According to Dr Wallace we update using previous not proposed step
567  // this way we do special proposal after adaptive after.
568  // This way we can shortcut and skip rest of proposal
569  if(!doSpecialStepProposal) return;
570 
572 }
void Randomize() _noexcept_
"Randomize" the parameters in the covariance class for the proposed step. Used the proposal kernel an...
void SpecialStepProposal()
Perform Special Step Proposal.
void CorrelateSteps() _noexcept_
Use Cholesky throw matrix for better step proposal.

◆ RandomConfiguration()

void ParameterHandlerBase::RandomConfiguration ( )

Randomly throw the parameters in their 1 sigma range.

Definition at line 504 of file ParameterHandlerBase.cpp.

504  {
505 // *************************************
506  // Have the 1 sigma for each parameter in each covariance class, sweet!
507  // Don't want to change the prior array because that's what determines our likelihood
508  // Want to change the _fPropVal, _fCurrVal, _fPreFitValue
509  // _fPreFitValue and the others will already be set
510  for (int i = 0; i < _fNumPar; ++i) {
511  // Check if parameter is fixed first: if so don't randomly throw
512  if (IsParameterFixed(i)) continue;
513  // Check that the sigma range is larger than the parameter range
514  // If not, throw in the valid parameter range instead
515  const double paramrange = _fUpBound[i] - _fLowBound[i];
516  const double sigma = sqrt((*covMatrix)(i,i));
517  double throwrange = sigma;
518  if (paramrange < sigma) throwrange = paramrange;
519 
520  _fPropVal[i] = _fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange;
521  // Try again if we the initial parameter proposal falls outside of the range of the parameter
522  int throws = 0;
523  while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
524  if (throws > 1000) {
525  MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
526  MACH3LOG_WARN("Matrix: {}", matrixName);
527  MACH3LOG_WARN("Param: {}", _fNames[i]);
528  throw MaCh3Exception(__FILE__ , __LINE__ );
529  }
530  _fPropVal[i] = _fPreFitValue[i] + random_number[0]->Gaus(0, 1)*throwrange;
531  throws++;
532  }
533  MACH3LOG_INFO("Setting current step in {} param {} = {} from {}", matrixName, i, _fPropVal[i], _fCurrVal[i]);
534  _fCurrVal[i] = _fPropVal[i];
535  }
536  if (pca) PCAObj->TransferToPCA();
537 }
std::vector< std::string > _fNames
ETA _fNames is set automatically in the covariance class to be something like xsec_i,...

◆ 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 598 of file ParameterHandlerBase.cpp.

598  {
599 // ************************************************
600  if (!pca) {
601  //KS: By multithreading here we gain at least factor 2 with 8 threads with ND only fit
602  #ifdef MULTITHREAD
603  #pragma omp parallel for
604  #endif
605  for (int i = 0; i < _fNumPar; ++i) {
606  // If parameter isn't fixed
607  if (!IsParameterFixed(i) > 0.0) {
608  randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
609  // If parameter IS fixed
610  } else {
611  randParams[i] = 0.0;
612  }
613  } // end for
614  // If we're in the PCA basis we instead throw parameters there (only _fNumParPCA parameter)
615  } else {
616  // Scale the random parameters by the sqrt of eigen values for the throw
617  #ifdef MULTITHREAD
618  #pragma omp parallel for
619  #endif
620  for (int i = 0; i < PCAObj->GetNumberPCAedParameters(); ++i)
621  {
622  // If parameter IS fixed or out of bounds
623  if (PCAObj->IsParameterFixedPCA(i)) {
624  randParams[i] = 0.0;
625  } else {
626  randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0,1);
627  }
628  }
629  }
630 }
int GetThreadIndex()
thread index inside parallel loop
Definition: Monitor.h:86

◆ ReserveMemory()

void ParameterHandlerBase::ReserveMemory ( const int  size)
protected

Initialise vectors with parameters information.

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

Definition at line 389 of file ParameterHandlerBase.cpp.

389  {
390 // ********************************************
391  _fNames = std::vector<std::string>(SizeVec);
392  _fFancyNames = std::vector<std::string>(SizeVec);
393  _fPreFitValue = std::vector<double>(SizeVec);
394  _fError = std::vector<double>(SizeVec);
395  _fCurrVal = std::vector<double>(SizeVec);
396  _fPropVal = std::vector<double>(SizeVec);
397  _fLowBound = std::vector<double>(SizeVec);
398  _fUpBound = std::vector<double>(SizeVec);
399  _fFlatPrior = std::vector<bool>(SizeVec);
400  _fIndivStepScale = std::vector<double>(SizeVec);
401  _fSampleNames = std::vector<std::vector<std::string>>(_fNumPar);
402 
403  corr_throw = new double[SizeVec];
404  // set random parameter vector (for correlated steps)
405  randParams = new double[SizeVec];
406 
407  // Set the defaults to true
408  for(int i = 0; i < SizeVec; i++) {
409  _fPreFitValue.at(i) = 1.;
410  _fError.at(i) = 1.;
411  _fCurrVal.at(i) = 0.;
412  _fPropVal.at(i) = 0.;
413  _fLowBound.at(i) = -999.99;
414  _fUpBound.at(i) = 999.99;
415  _fFlatPrior.at(i) = false;
416  _fIndivStepScale.at(i) = 1.;
417  corr_throw[i] = 0.0;
418  randParams[i] = 0.0;
419  }
420 
421  _fGlobalStepScale = 1.0;
422 }

◆ ResetIndivStepScale()

void ParameterHandlerBase::ResetIndivStepScale ( )

Adaptive Step Tuning Stuff.

Definition at line 1007 of file ParameterHandlerBase.cpp.

1007  {
1008 // ********************************************
1009  std::vector<double> stepScales(_fNumPar);
1010  for (int i = 0; i <_fNumPar; i++) {
1011  stepScales[i] = 1.;
1012  }
1013  _fGlobalStepScale = 1.0;
1014  SetIndivStepScale(stepScales);
1015 }
void SetIndivStepScale(const int ParameterIndex, const double StepScale)
DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constr...

◆ RetPointer()

const double* ParameterHandlerBase::RetPointer ( const int  iParam)
inline

DB Pointer return to param position.

Parameters
iParamThe 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 233 of file ParameterHandlerBase.h.

233 {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 193 of file ParameterHandlerBase.h.

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

◆ SaveUpdatedMatrixConfig()

void ParameterHandlerBase::SaveUpdatedMatrixConfig ( )

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

Definition at line 1210 of file ParameterHandlerBase.cpp.

1210  {
1211 // ********************************************
1212  if (!_fYAMLDoc)
1213  {
1214  MACH3LOG_CRITICAL("Yaml node hasn't been initialised for matrix {}, something is not right", matrixName);
1215  MACH3LOG_CRITICAL("I am not throwing error but should be investigated");
1216  return;
1217  }
1218 
1219  YAML::Node copyNode = _fYAMLDoc;
1220  int i = 0;
1221 
1222  for (YAML::Node param : copyNode["Systematics"])
1223  {
1224  //KS: Feel free to update it, if you need updated prefit value etc
1225  param["Systematic"]["StepScale"]["MCMC"] = MaCh3Utils::FormatDouble(_fIndivStepScale[i], 4);
1226  i++;
1227  }
1228  // Save the modified node to a file
1229  std::ofstream fout("Modified_Matrix.yaml");
1230  fout << copyNode;
1231  fout.close();
1232 }
std::string FormatDouble(const double value, const int precision)
Convert double into string for precision, useful for playing with yaml if you don't want to have in c...

◆ 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_namename of file matrix lives in
matrix_namename of matrix in file
means_namename of means vec in file

◆ SpecialStepProposal()

void ParameterHandlerBase::SpecialStepProposal ( )
protected

Perform Special Step Proposal.

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

Definition at line 575 of file ParameterHandlerBase.cpp.

575  {
576 // ************************************************
578 
579  // HW It should now automatically set dcp to be with [-pi, pi]
580  for (size_t i = 0; i < CircularBoundsIndex.size(); ++i) {
581  const int index = CircularBoundsIndex[i];
582  if(!IsParameterFixed(index))
583  CircularParBounds(index, CircularBoundsValues[i].first, CircularBoundsValues[i].second);
584  }
585 
586  // Okay now we've done the standard steps, we can add in our nice flips hierarchy flip first
587  for (size_t i = 0; i < FlipParameterIndex.size(); ++i) {
588  const int index = FlipParameterIndex[i];
589  if(!IsParameterFixed(index))
591  }
592 }
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.
void FlipParameterValue(const int index, const double FlipPoint)
KS: Flip parameter around given value, for example mass ordering around 0.

◆ 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 451 of file ParameterHandlerBase.cpp.

451  {
452 // *************************************
453  // First draw new randParams
454  Randomize();
455 
457 
458  // KS: We use PCA very rarely on top PCA functionality isn't implemented for this function.
459  // Use __builtin_expect to give compiler a hint which option is more likely, which should help
460  // with better optimisation. This isn't critical but more to have example
461  if (__builtin_expect(!pca, 1)) {
462  #ifdef MULTITHREAD
463  #pragma omp parallel for
464  #endif
465  for (int i = 0; i < _fNumPar; ++i) {
466  // Check if parameter is fixed first: if so don't randomly throw
467  if (IsParameterFixed(i)) continue;
468 
469  _fPropVal[i] = _fPreFitValue[i] + corr_throw[i];
470  int throws = 0;
471  // Try again if we the initial parameter proposal falls outside of the range of the parameter
472  while (_fPropVal[i] > _fUpBound[i] || _fPropVal[i] < _fLowBound[i]) {
473  randParams[i] = random_number[M3::GetThreadIndex()]->Gaus(0, 1);
474  const double corr_throw_single = M3::MatrixVectorMultiSingle(throwMatrixCholDecomp, randParams, _fNumPar, i);
475  _fPropVal[i] = _fPreFitValue[i] + corr_throw_single;
476  if (throws > 10000)
477  {
478  //KS: Since we are multithreading there is danger that those messages
479  //will be all over the place, small price to pay for faster code
480  MACH3LOG_WARN("Tried {} times to throw parameter {} but failed", throws, i);
481  MACH3LOG_WARN("Matrix: {}", matrixName);
482  MACH3LOG_WARN("Param: {}", _fNames[i]);
483  MACH3LOG_WARN("Setting _fPropVal: {} to {}", _fPropVal[i], _fPreFitValue[i]);
484  MACH3LOG_WARN("I live at {}:{}", __FILE__, __LINE__);
485  _fPropVal[i] = _fPreFitValue[i];
486  //throw MaCh3Exception(__FILE__ , __LINE__ );
487  }
488  throws++;
489  }
490  _fCurrVal[i] = _fPropVal[i];
491  }
492  }
493  else
494  {
495  PCAObj->ThrowParameters(random_number, throwMatrixCholDecomp,
498  } // end if pca
499 }
double MatrixVectorMultiSingle(double **_restrict_ matrix, const double *_restrict_ vector, const int Length, const int i)
KS: Custom function to perform multiplication of matrix and single element which is thread safe.

◆ ThrowParCurr()

void ParameterHandlerBase::ThrowParCurr ( const double  mag = 1.)

Helper function to throw the current parameter by mag sigma. Can study bias in MCMC with this; put different starting parameters.

Definition at line 711 of file ParameterHandlerBase.cpp.

711  {
712 // ********************************************
713  Randomize();
714  if (!pca) {
715  // Get the correlated throw vector
717  // The number of sigmas to throw
718  // Should probably have this as a default parameter input to the function instead
719  for (int i = 0; i < _fNumPar; i++) {
720  if (!IsParameterFixed(i) > 0.){
721  _fCurrVal[i] = corr_throw[i]*mag;
722  }
723  }
724  } else {
725  PCAObj->ThrowParCurr(mag, randParams);
726  }
727 }

◆ ThrowParProp()

void ParameterHandlerBase::ThrowParProp ( const double  mag = 1.)

Throw the proposed parameter by mag sigma. Should really just have the user specify this throw by having argument double.

Definition at line 693 of file ParameterHandlerBase.cpp.

693  {
694 // ********************************************
695  Randomize();
696  if (!pca) {
697  // Make the correlated throw
699  // Number of sigmas we throw
700  for (int i = 0; i < _fNumPar; i++) {
701  if (!IsParameterFixed(i) > 0.)
702  _fPropVal[i] = _fCurrVal[i] + corr_throw[i]*mag;
703  }
704  } else {
705  PCAObj->ThrowParProp(mag, randParams);
706  }
707 }

◆ ToggleFixAllParameters()

void ParameterHandlerBase::ToggleFixAllParameters ( )

fix parameters at prior values

Definition at line 896 of file ParameterHandlerBase.cpp.

896  {
897 // ********************************************
898  // fix or unfix all parameters by multiplying by -1
899  if(!pca) {
900  for (int i = 0; i < _fNumPar; i++) _fError[i] *= -1.0;
901  } else{
902  PCAObj->ToggleFixAllParameters();
903  }
904 }

◆ ToggleFixParameter() [1/2]

void ParameterHandlerBase::ToggleFixParameter ( const int  i)

fix parameter at prior values

Parameters
iParameter index

Definition at line 907 of file ParameterHandlerBase.cpp.

907  {
908 // ********************************************
909  if(!pca) {
910  if (i > _fNumPar) {
911  MACH3LOG_ERROR("Can't {} for parameter {} because size of covariance ={}", __func__, i, _fNumPar);
912  MACH3LOG_ERROR("Fix this in your config file please!");
913  throw MaCh3Exception(__FILE__ , __LINE__ );
914  } else {
915  _fError[i] *= -1.0;
916  MACH3LOG_INFO("Setting {}(parameter {}) to fixed at {}", GetParFancyName(i), i, _fCurrVal[i]);
917  }
918  } else {
919  PCAObj->ToggleFixParameter(i, _fNames);
920  }
921 }

◆ ToggleFixParameter() [2/2]

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

Fix parameter at prior values.

Parameters
nameName of parameter you want to fix

Definition at line 924 of file ParameterHandlerBase.cpp.

924  {
925 // ********************************************
926  const int Index = GetParIndex(name);
927  if(Index != M3::_BAD_INT_) {
928  ToggleFixParameter(Index);
929  return;
930  }
931 
932  MACH3LOG_WARN("I couldn't find parameter with name {}, therefore will not fix it", name);
933 }

◆ UpdateAdaptiveCovariance()

void ParameterHandlerBase::UpdateAdaptiveCovariance ( )

Method to update adaptive MCMC [12].

Definition at line 1102 of file ParameterHandlerBase.cpp.

1102  {
1103 // ********************************************
1104  // Updates adaptive matrix
1105  // First we update the total means
1106 
1107  // Skip this if we're at a large number of steps
1108  if(AdaptiveHandler->SkipAdaption()) {
1109  AdaptiveHandler->IncrementNSteps();
1110  return;
1111  }
1112 
1113  // Call main adaption function
1114  AdaptiveHandler->UpdateAdaptiveCovariance();
1115 
1116  //This is likely going to be the slow bit!
1117  if(AdaptiveHandler->IndivStepScaleAdapt()) {
1119  }
1120 
1121  if(AdaptiveHandler->UpdateMatrixAdapt()) {
1122  TMatrixDSym* update_matrix = static_cast<TMatrixDSym*>(AdaptiveHandler->GetAdaptiveCovariance()->Clone());
1123  UpdateThrowMatrix(update_matrix); //Now we update and continue!
1124  //Also Save the adaptive to file
1125  AdaptiveHandler->SaveAdaptiveToFile(AdaptiveHandler->GetOutFileName(), GetName());
1126  }
1127 
1128  AdaptiveHandler->IncrementNSteps();
1129 }
void UpdateThrowMatrix(TMatrixDSym *cov)

◆ UpdateThrowMatrix()

void ParameterHandlerBase::UpdateThrowMatrix ( TMatrixDSym *  cov)

Definition at line 1053 of file ParameterHandlerBase.cpp.

1053  {
1054 // ********************************************
1055  delete throwMatrix;
1056  throwMatrix = nullptr;
1057  SetThrowMatrix(cov);
1058 }

Member Data Documentation

◆ _fCurrVal

std::vector<double> ParameterHandlerBase::_fCurrVal
protected

Current value of the parameter.

Definition at line 459 of file ParameterHandlerBase.h.

◆ _fError

std::vector<double> ParameterHandlerBase::_fError
protected

Prior error on the parameter.

Definition at line 463 of file ParameterHandlerBase.h.

◆ _fFancyNames

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

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

Definition at line 451 of file ParameterHandlerBase.h.

◆ _fFlatPrior

std::vector<bool> ParameterHandlerBase::_fFlatPrior
protected

Whether to apply flat prior or not.

Definition at line 471 of file ParameterHandlerBase.h.

◆ _fGlobalStepScale

double ParameterHandlerBase::_fGlobalStepScale
protected

Global step scale applied to all params in this class.

Definition at line 443 of file ParameterHandlerBase.h.

◆ _fIndivStepScale

std::vector<double> ParameterHandlerBase::_fIndivStepScale
protected

Individual step scale used by MCMC algorithm.

Definition at line 469 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 465 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 xsec_i, this is currently to make things compatible with the Diagnostic tools.

Definition at line 449 of file ParameterHandlerBase.h.

◆ _fNumPar

int ParameterHandlerBase::_fNumPar
protected

Number of systematic parameters.

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

◆ _fPropVal

std::vector<double> ParameterHandlerBase::_fPropVal
protected

Proposed value of the parameter.

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

◆ _fYAMLDoc

YAML::Node ParameterHandlerBase::_fYAMLDoc
protected

Stores config describing systematics.

Definition at line 453 of file ParameterHandlerBase.h.

◆ AdaptiveHandler

std::unique_ptr<adaptive_mcmc::AdaptiveMCMCHandler> ParameterHandlerBase::AdaptiveHandler
protected

Struct containing information about adaption.

Definition at line 488 of file ParameterHandlerBase.h.

◆ CircularBoundsIndex

std::vector<int> ParameterHandlerBase::CircularBoundsIndex
protected

Indices of parameters with circular bounds.

Definition at line 497 of file ParameterHandlerBase.h.

◆ CircularBoundsValues

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

Circular bounds for each parameter (lower, upper)

Definition at line 499 of file ParameterHandlerBase.h.

◆ corr_throw

double* ParameterHandlerBase::corr_throw
protected

Result of multiplication of Cholesky matrix and randParams.

Definition at line 441 of file ParameterHandlerBase.h.

◆ covMatrix

TMatrixDSym* ParameterHandlerBase::covMatrix
protected

The covariance matrix.

Definition at line 429 of file ParameterHandlerBase.h.

◆ doSpecialStepProposal

bool ParameterHandlerBase::doSpecialStepProposal
protected

Check if any of special step proposal were enabled.

Definition at line 421 of file ParameterHandlerBase.h.

◆ FlipParameterIndex

std::vector<int> ParameterHandlerBase::FlipParameterIndex
protected

Indices of parameters with flip symmetry.

Definition at line 493 of file ParameterHandlerBase.h.

◆ FlipParameterPoint

std::vector<double> ParameterHandlerBase::FlipParameterPoint
protected

Central points around which parameters are flipped.

Definition at line 495 of file ParameterHandlerBase.h.

◆ inputFile

const std::string ParameterHandlerBase::inputFile
protected

The input root file we read in.

Definition at line 424 of file ParameterHandlerBase.h.

◆ invCovMatrix

TMatrixDSym* ParameterHandlerBase::invCovMatrix
protected

The inverse covariance matrix.

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

◆ matrixName

std::string ParameterHandlerBase::matrixName
protected

Name of cov matrix.

Definition at line 427 of file ParameterHandlerBase.h.

◆ pca

bool ParameterHandlerBase::pca
protected

perform PCA or not

Definition at line 481 of file ParameterHandlerBase.h.

◆ PCAObj

std::unique_ptr<PCAHandler> ParameterHandlerBase::PCAObj
protected

Struct containing information about PCA.

Definition at line 486 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 446 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 436 of file ParameterHandlerBase.h.

◆ randParams

double* ParameterHandlerBase::randParams
protected

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

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

◆ Tunes

std::unique_ptr<ParameterTunes> ParameterHandlerBase::Tunes
protected

Struct containing information about adaption.

Definition at line 490 of file ParameterHandlerBase.h.

◆ use_adaptive

bool ParameterHandlerBase::use_adaptive
protected

Are we using AMCMC?

Definition at line 483 of file ParameterHandlerBase.h.


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