tidy: Step Scale Sanity chekc et. al.

.github/workflows/CIValidations.yml

@@ -25,7 +25,7 @@ jobs:
             test_1: ./CIValidations/SplineValidations
             test_2: ./CIValidations/SamplePDFValidations
             test_3: ./CIValidations/NuOscillatorInterfaceValidations
-            test_4: empty
+            test_4: ./CIValidations/LLHValidation
             test_5: empty
             test_6: empty
             test_7: empty

Doc/bibliography.bib

@@ -223,3 +223,26 @@ @manual{Hartig2024WAIC
   url          = {https://search.r-project.org/CRAN/refmans/BayesianTools/html/WAIC.html},
   institution  = {CRAN},
 }
+
+@article{luengo2020survey,
+  author = {Luengo, D. and Martino, L. and Bugallo, M. and others},
+  title = {A survey of Monte Carlo methods for parameter estimation},
+  journal = {EURASIP Journal on Advances in Signal Processing},
+  year = {2020},
+  doi = {10.1186/s13634-020-00675-6},
+  url = {https://doi.org/10.1186/s13634-020-00675-6},
+  note = {Relevant discussion in Section 3.2.1},
+  published = {29 May 2020}
+}
+
+@article{haario2001adaptive,
+  author       = {Heikki Haario and Eero Saksman and Johanna Tamminen},
+  title        = {An adaptive Metropolis algorithm},
+  journal      = {Bernoulli},
+  volume       = {7},
+  number       = {2},
+  pages        = {223--242},
+  year         = {2001},
+  month        = {April},
+  doi          = {10.3150/bj/1080222083}
+}

covariance/AdaptiveMCMCHandler.h

@@ -7,7 +7,8 @@
 namespace adaptive_mcmc{
 
 /// @brief Contains information about adaptive covariance matrix
-/// @see An adaptive Metropolis algorithm, H.Haario et al., 2001 for more info!
+/// @cite haario2001adaptive
+/// @author Henry Wallace
 /// @details struct encapsulating all adaptive MCMC information
 class AdaptiveMCMCHandler{
  public:
@@ -16,7 +17,7 @@ class AdaptiveMCMCHandler{
   AdaptiveMCMCHandler();
 
   /// @brief Destructor
-  ~AdaptiveMCMCHandler();
+  virtual ~AdaptiveMCMCHandler();
 
   /// @brief Print all class members
   void Print();
@@ -46,7 +47,7 @@ class AdaptiveMCMCHandler{
                               const int Npars);
 
   /// @brief Method to update adaptive MCMC
-  /// @see https://projecteuclid.org/journals/bernoulli/volume-7/issue-2/An-adaptive-Metropolis-algorithm/bj/1080222083.full
+  /// @cite haario2001adaptive
   /// @param _fCurrVal Value of each parameter necessary for updating throw matrix
   void UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int Npars);
 

covariance/PCAHandler.h

@@ -33,7 +33,7 @@ class PCAHandler{
   PCAHandler();
 
   /// @brief Destructor
-  ~PCAHandler();
+  virtual ~PCAHandler();
 
   /// @brief CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold
   void ConstructPCA(TMatrixDSym * covMatrix, const int firstPCAd, const int lastPCAd, const double eigen_thresh, int& _fNumParPCA);

covariance/covarianceBase.cpp

@@ -60,16 +60,11 @@ covarianceBase::~covarianceBase(){
   if (covMatrix != nullptr) delete covMatrix;
   if (invCovMatrix != nullptr) delete invCovMatrix;
   if (throwMatrix_CholDecomp != nullptr) delete throwMatrix_CholDecomp;
-
-  for(int i = 0; i < _fNumPar; i++)
-  {
-    delete[] InvertCovMatrix[i];
+  if (throwMatrix != nullptr) delete throwMatrix;
+  for(int i = 0; i < _fNumPar; i++) {
     delete[] throwMatrixCholDecomp[i];
   }
-  delete[] InvertCovMatrix;
   delete[] throwMatrixCholDecomp;
-
-  if (throwMatrix != nullptr) delete throwMatrix;
 }
 
 // ********************************************
@@ -112,7 +107,6 @@ void covarianceBase::ConstructPCA() {
 // ********************************************
 void covarianceBase::init(std::string name, std::string file) {
 // ********************************************
-
   // Set the covariance matrix from input ROOT file (e.g. flux, ND280, NIWG)
   TFile *infile = new TFile(file.c_str(), "READ");
   if (infile->IsZombie()) {
@@ -143,20 +137,15 @@ void covarianceBase::init(std::string name, std::string file) {
   // Set the covariance matrix
   _fNumPar = CovMat->GetNrows();
 
-  InvertCovMatrix = new double*[_fNumPar]();
+  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++)
-  {
-    InvertCovMatrix[i] = new double[_fNumPar]();
+  for(int i = 0; i < _fNumPar; i++) {
     throwMatrixCholDecomp[i] = new double[_fNumPar]();
-    for (int j = 0; j < _fNumPar; j++)
-    {
-      InvertCovMatrix[i][j] = 0.;
+    for (int j = 0; j < _fNumPar; j++) {
       throwMatrixCholDecomp[i][j] = 0.;
     }
   }
-
   setName(name);
   MakePosDef(CovMat);
   setCovMatrix(CovMat);
@@ -172,7 +161,6 @@ void covarianceBase::init(std::string name, std::string file) {
 
   MACH3LOG_INFO("Created covariance matrix named: {}", getName());
   MACH3LOG_INFO("from file: {}", file);
-
   delete infile;
 }
 
@@ -206,25 +194,17 @@ void covarianceBase::init(const std::vector<std::string>& YAMLFile) {
 
   use_adaptive = false;
 
-  InvertCovMatrix = new double*[_fNumPar]();
+  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++)
-  {
-    InvertCovMatrix[i] = new double[_fNumPar]();
+  for(int i = 0; i < _fNumPar; i++) {
     throwMatrixCholDecomp[i] = new double[_fNumPar]();
-    for (int j = 0; j < _fNumPar; j++)
-    {
-      InvertCovMatrix[i][j] = 0.;
+    for (int j = 0; j < _fNumPar; j++) {
       throwMatrixCholDecomp[i][j] = 0.;
     }
   }
-
   ReserveMemory(_fNumPar);
 
   TMatrixDSym* _fCovMatrix = new TMatrixDSym(_fNumPar);
-  //_fDetString = std::vector<std::string>(_fNumPar);
-
   int i = 0;
   std::vector<std::map<std::string,double>> Correlations(_fNumPar);
   std::map<std::string, int> CorrNamesMap;
@@ -239,7 +219,10 @@ void covarianceBase::init(const std::vector<std::string>& YAMLFile) {
     _fGenerated[i] = (param["Systematic"]["ParameterValues"]["Generated"].as<double>());
     _fIndivStepScale[i] = (param["Systematic"]["StepScale"]["MCMC"].as<double>());
     _fError[i] = (param["Systematic"]["Error"].as<double>());
-
+    if(_fError[i] <= 0) {
+      MACH3LOG_ERROR("Error for param {}({}) is negative and eqaul to {}", _fFancyNames[i], i, _fError[i]);
+      throw MaCh3Exception(__FILE__ , __LINE__ );
+    }
     //ETA - a bit of a fudge but works
     std::vector<double> TempBoundsVec = param["Systematic"]["ParameterBounds"].as<std::vector<double>>();
     _fLowBound[i] = TempBoundsVec[0];
@@ -265,12 +248,14 @@ void covarianceBase::init(const std::vector<std::string>& YAMLFile) {
     }
     i++;
   }
-
-  //ETA
-  //Now that we've been through all systematic let's fill the covmatrix
+  if(i != _fNumPar) {
+    MACH3LOG_CRITICAL("Inconsitnent 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];
+  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;
@@ -367,7 +352,7 @@ void covarianceBase::ReserveMemory(const int SizeVec) {
   for(int i = 0; i < SizeVec; i++) {
     _fGenerated.at(i) = 1.;
     _fPreFitValue.at(i) = 1.;
-    _fError.at(i) = 0.;
+    _fError.at(i) = 1.;
     _fCurrVal.at(i) = 0.;
     _fPropVal.at(i) = 0.;
     _fLowBound.at(i) = -999.99;
@@ -820,8 +805,7 @@ double covarianceBase::GetLikelihood() {
   // Default behaviour is to reject negative values + do std llh calculation
   const int NOutside = CheckBounds();
 
-  if(NOutside > 0)
-    return NOutside*_LARGE_LOGL_;
+  if(NOutside > 0) return NOutside*_LARGE_LOGL_;
 
   return CalcLikelihood();
 }
@@ -904,21 +888,24 @@ void covarianceBase::SetBranches(TTree &tree, bool SaveProposal) {
 // ********************************************
 void covarianceBase::setStepScale(const double scale) {
 // ********************************************
-  if(scale == 0)
-  {
-    MACH3LOG_ERROR("You are trying so set StepScale to 0 this will not work");
+  if(scale <= 0) {
+    MACH3LOG_ERROR("You are trying so set StepScale to 0 or negative this will not work");
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   MACH3LOG_INFO("{} setStepScale() = {}", getName(), scale);
+  const double SuggestedScale = 2.38*2.38/_fNumPar;
+  if(std::fabs(scale - SuggestedScale)/SuggestedScale > 1) {
+    MACH3LOG_WARN("Defined Global StepScale is {}, while suggested suggested {}", scale, SuggestedScale);
+  }
+
   _fGlobalStepScale = scale;
 }
 
 // ********************************************
 void covarianceBase::toggleFixAllParameters() {
 // ********************************************
   // fix or unfix all parameters by multiplying by -1
-  if(!pca)
-  {
+  if(!pca) {
     for (int i = 0; i < _fNumPar; i++) _fError[i] *= -1.0;
   } else{
      for (int i = 0; i < _fNumParPCA; i++) fParSigma_PCA[i] *= -1.0;
@@ -1171,28 +1158,25 @@ void covarianceBase::updateThrowMatrix(TMatrixDSym *cov){
 // HW : Here be adaption
 void covarianceBase::initialiseAdaption(const YAML::Node& adapt_manager){
 // ********************************************
-  // need to cast matrixName to string
-  std::string matrix_name_str(matrixName);
-
   // Now we read the general settings [these SHOULD be common across all matrices!]
-  bool success = AdaptiveHandler.InitFromConfig(adapt_manager, matrix_name_str, getNpars());
+  bool success = AdaptiveHandler.InitFromConfig(adapt_manager, matrixName, getNpars());
   if(!success) return;
   AdaptiveHandler.Print();
 
   // Next let"s check for external matrices
   // We"re going to grab this info from the YAML manager
-  if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["UseExternalMatrix"], false)) {
-    MACH3LOG_WARN("Not using external matrix for {}, initialising adaption from scratch", matrix_name_str);
+  if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["UseExternalMatrix"], false)) {
+    MACH3LOG_WARN("Not using external matrix for {}, initialising adaption from scratch", matrixName);
     // If we don't have a covariance matrix to start from for adaptive tune we need to make one!
     use_adaptive = true;
     AdaptiveHandler.CreateNewAdaptiveCovariance(_fNumPar);
     return;
   }
 
   // Finally, we accept that we want to read the matrix from a file!
-  std::string external_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMatrixFileName"], "");
-  std::string external_matrix_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMatrixName"], "");
-  std::string external_mean_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMeansName"], "");
+  std::string external_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixFileName"], "");
+  std::string external_matrix_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMatrixName"], "");
+  std::string external_mean_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrixName]["ExternalMeansName"], "");
 
   AdaptiveHandler.SetThrowMatrixFromFile(external_file_name, external_matrix_name, external_mean_name, use_adaptive, _fNumPar);
   setThrowMatrix(AdaptiveHandler.adaptive_covariance);

covariance/covarianceBase.h

@@ -69,6 +69,7 @@ class covarianceBase {
   void SetBranches(TTree &tree, bool SaveProposal = false);
   /// @brief Set global step scale for covariance object
   /// @param scale Value of global step scale
+  /// @cite luengo2020survey
   void setStepScale(const double scale);
   /// @brief DB Function to set fIndivStepScale from a vector (Can be used from execs and inside covariance constructors)
   /// @param ParameterIndex Parameter Index
@@ -389,7 +390,7 @@ class covarianceBase {
   /// The inverse covariance matrix
   TMatrixDSym *invCovMatrix;
   /// KS: Same as above but much faster as TMatrixDSym cache miss
-  double **InvertCovMatrix;
+  std::vector<std::vector<double>> InvertCovMatrix;
 
   /// KS: Set Random numbers for each thread so each thread has different seed
   std::vector<std::unique_ptr<TRandom3>> random_number;
@@ -461,7 +462,7 @@ class covarianceBase {
   /// Matrix which we use for step proposal after Cholesky decomposition
   TMatrixD* throwMatrix_CholDecomp;
   /// Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss
-  double **throwMatrixCholDecomp;
+  double** throwMatrixCholDecomp;
 
   /// Are we using AMCMC?
   bool use_adaptive;

mcmc/FitterBase.cpp

@@ -956,7 +956,6 @@ void FitterBase::Run2DLLHScan() {
           MACH3LOG_INFO("lower {} = {:.2f}", i, lower_y);
           MACH3LOG_INFO("upper {} = {:.2f}", i, upper_y);
           MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
-
         }
 
         // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables