Merge pull request #80 from mach3-software/feature_AapationTidy

Adaptation Handler

.mailmap

@@ -12,3 +12,5 @@ Kamil Skwarczynski <skwarczynskikamil@gmail.com> Kamil Skwarczynski <kamilskw@cd
 # Henry Wallace
 Henry Wallace <henryisrael08@gmail.com> henry-israel <henryisrael08@gmail.com>
 Henry Wallace <henryisrael08@gmail.com> Henry Wallace <67589487+henry-wallace-phys@users.noreply.github.com>
+Henry Wallace <henryisrael08@gmail.com> Henry Wallace <hwallace@linappserv6.pp.rhul.ac.uk>
+Henry Wallace <henryisrael08@gmail.com> Henry Wallace <henry.wallace@rhul.ac.uk>

CMakeLists.txt

@@ -314,3 +314,6 @@ if(NOT TARGET uninstall)
   add_custom_target(uninstall
     COMMAND ${CMAKE_COMMAND} -P ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake)
 endif()
+
+# KS: Configure the Doxygen input file, this is to ensure whenever we update MaCh3 version Doxyfile will have same version.
+configure_file(${CMAKE_SOURCE_DIR}/cmake/Templates/Doxyfile.in ${CMAKE_SOURCE_DIR}/Doc/Doxyfile @ONLY)

SECURITY.md

@@ -0,0 +1,13 @@
+# Security Policy
+
+## Supported Versions
+
+Security updates are applied only to the latest release.
+
+## Reporting a Vulnerability
+
+If you have discovered a security vulnerability in this project, please report it privately. **Do not disclose it as a public issue.** This gives us time to work with you to fix the issue before public exposure, reducing the chance that the exploit will be used before a patch is released.
+
+Please disclose it at [security advisory](https://github.com/mach3-software/mach3/security/advisories/new).
+
+This project is maintained by a team of volunteers on a reasonable-effort basis. As such, vulnerabilities will be disclosed in a best effort base.

covariance/AdaptiveMCMCHandler.cpp

@@ -0,0 +1,231 @@
+#include "covariance/AdaptiveMCMCHandler.h"
+
+namespace adaptive_mcmc{
+
+// ********************************************
+AdaptiveMCMCHandler::AdaptiveMCMCHandler() {
+// ********************************************
+  start_adaptive_throw  = 0;
+  start_adaptive_update = 0;
+  end_adaptive_update   = 1;
+  adaptive_update_step  = 1000;
+
+  par_means = {};
+  adaptive_covariance = nullptr;
+}
+
+// ********************************************
+AdaptiveMCMCHandler::~AdaptiveMCMCHandler() {
+// ********************************************
+  if(adaptive_covariance != nullptr) {
+    delete adaptive_covariance;
+  }
+}
+
+
+// ********************************************
+void AdaptiveMCMCHandler::InitFromConfig(const YAML::Node& adapt_manager, const std::string& matrix_name_str, const int Npars) {
+// ********************************************
+  //  setAdaptionDefaults();
+  if(GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str], "")==""){
+    MACH3LOG_WARN("Adaptive Settings not found for {}, this is fine if you don't want adaptive MCMC", matrix_name_str);
+    return;
+  }
+
+  // We"re going to grab this info from the YAML manager
+  if(GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["DoAdaption"], false)) {
+    MACH3LOG_WARN("Not using adaption for {}", matrix_name_str);
+    return;
+  }
+
+  start_adaptive_throw  = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionStartThrow"], 10);
+  start_adaptive_update = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionStartUpdate"], 0);
+  end_adaptive_update   = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionEndUpdate"], 10000);
+  adaptive_update_step  = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionUpdateStep"], 100);
+
+
+  // We also want to check for "blocks" by default all parameters "know" about each other
+  // but we can split the matrix into independent block matrices
+
+  // We"ll set a dummy variable here
+  auto matrix_blocks = GetFromManager<std::vector<std::vector<int>>>(adapt_manager["AdaptionOptions"]["Settings"][matrix_name_str]["AdaptionUpdateStep"], {{}});
+
+  SetAdaptiveBlocks(matrix_blocks, Npars);
+}
+
+// ********************************************
+void AdaptiveMCMCHandler::CreateNewAdaptiveCovariance(const int Npars) {
+// ********************************************
+  adaptive_covariance = new TMatrixDSym(Npars);
+  adaptive_covariance->Zero();
+  par_means = std::vector<double>(Npars, 0);
+}
+
+
+// ********************************************
+void AdaptiveMCMCHandler::SetAdaptiveBlocks(std::vector<std::vector<int>> block_indices, const int Npars) {
+// ********************************************
+  /*
+   *   In order to adapt efficient we want to setup our throw matrix to be a serious of block-diagonal (ish) matrices
+   *
+   *   To do this we set sub-block in the config by parameter index. For example having
+   *   [[0,4],[4, 6]] in your config will set up two blocks one with all indices 0<=i<4 and the other with 4<=i<6
+   */
+  // Set up block regions
+  adapt_block_matrix_indices = std::vector<int>(Npars, 0);
+
+  // Should also make a matrix of block sizes
+  adapt_block_sizes = std::vector<int>((int)block_indices.size()+1, 0);
+  adapt_block_sizes[0] = Npars;
+
+  if(block_indices.size()==0 || block_indices[0].size()==0) return;
+
+  // Now we loop over our blocks
+  for(int iblock=0; iblock<(int)block_indices.size(); iblock++){
+    // Loop over blocks in the block
+    for(int isubblock=0; isubblock<(int)block_indices[iblock].size()-1; isubblock+=2){
+      int block_lb = block_indices[iblock][isubblock];
+      int block_ub = block_indices[iblock][isubblock+1];
+
+      //std::cout<<block_lb<<" "<<block_ub<<std::endl;
+      if(block_lb > Npars || block_ub > Npars){
+        MACH3LOG_ERROR("Cannot set matrix block with edges {}, {} for matrix of size {}",
+                       block_lb, block_ub, Npars);
+        throw MaCh3Exception(__FILE__, __LINE__);;
+      }
+      for(int ipar = block_lb; ipar < block_ub; ipar++){
+        adapt_block_matrix_indices[ipar] = iblock+1;
+        adapt_block_sizes[iblock+1] += 1;
+        adapt_block_sizes[0] -= 1;
+      }
+    }
+  }
+}
+
+
+// ********************************************
+//HW: Truly adaptive MCMC!
+void AdaptiveMCMCHandler::SaveAdaptiveToFile(const TString& outFileName, const TString& systematicName){
+// ********************************************
+  TFile* outFile = new TFile(outFileName, "UPDATE");
+  if(outFile->IsZombie()){
+    MACH3LOG_ERROR("Couldn't find {}", outFileName);
+    throw;
+  }
+  TVectorD* outMeanVec = new TVectorD((int)par_means.size());
+  for(int i = 0; i < (int)par_means.size(); i++){
+    (*outMeanVec)(i) = par_means[i];
+  }
+  outFile->cd();
+  adaptive_covariance->Write(systematicName+"_postfit_matrix");
+  outMeanVec->Write(systematicName+"_mean_vec");
+  outFile->Close();
+  delete outMeanVec;
+  delete outFile;
+}
+
+// ********************************************
+// HW : I would like this to be less painful to use!
+// First things first we need setters
+void AdaptiveMCMCHandler::SetThrowMatrixFromFile(const std::string& matrix_file_name,
+                                                 const std::string& matrix_name,
+                                                 const std::string& means_name,
+                                                 bool& use_adaptive,
+                                                 const int Npars) {
+// ********************************************
+  // Lets you set the throw matrix externally
+  // Open file
+  std::unique_ptr<TFile>matrix_file(new TFile(matrix_file_name.c_str()));
+  use_adaptive = true;
+
+  if(matrix_file->IsZombie()){
+    MACH3LOG_ERROR("Couldn't find {}", matrix_file_name);
+    throw;
+  }
+
+  // Next we grab our matrix
+  adaptive_covariance = static_cast<TMatrixDSym*>(matrix_file->Get(matrix_name.c_str()));
+  if(!adaptive_covariance){
+    MACH3LOG_ERROR("Couldn't find {} in {}", matrix_name, matrix_file_name);
+    throw;
+  }
+
+  // Finally we grab the means vector
+  TVectorD* means_vector = static_cast<TVectorD*>(matrix_file->Get(means_name.c_str()));
+
+  // This is fine to not exist!
+  if(means_vector){
+    // Yay our vector exists! Let's loop and fill it
+    // Should check this is done
+    if(means_vector->GetNrows()){
+      MACH3LOG_ERROR("External means vec size ({}) != matrix size ({})", means_vector->GetNrows(), Npars);
+      throw MaCh3Exception(__FILE__, __LINE__);
+    }
+
+    par_means = std::vector<double>(Npars);
+    for(int i = 0; i < Npars; i++){
+      par_means[i] = (*means_vector)(i);
+    }
+    MACH3LOG_INFO("Found Means in External File, Will be able to adapt");
+  }
+  // Totally fine if it doesn't exist, we just can't do adaption
+  else{
+    // We don't need a means vector, set the adaption=false
+    MACH3LOG_WARN("Cannot find means vector in {}, therefore I will not be able to adapt!", matrix_file_name);
+    use_adaptive = false;
+  }
+
+  matrix_file->Close();
+  MACH3LOG_INFO("Set up matrix from external file");
+}
+
+
+// ********************************************
+void AdaptiveMCMCHandler::UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int steps_post_burn, const int Npars) {
+// ********************************************
+  std::vector<double> par_means_prev = par_means;
+
+  #ifdef MULTITHREAD
+  #pragma omp parallel for
+  #endif
+  for(int iRow = 0; iRow < Npars; iRow++) {
+    par_means[iRow] = (_fCurrVal[iRow]+par_means[iRow]*steps_post_burn)/(steps_post_burn+1);
+  }
+
+  //Now we update the covariances using cov(x,y)=E(xy)-E(x)E(y)
+  #ifdef MULTITHREAD
+  #pragma omp parallel for
+  #endif
+  for(int irow = 0; irow < Npars; irow++){
+    int block = adapt_block_matrix_indices[irow];
+    // int scale_factor = 5.76/double(adapt_block_sizes[block]);
+    for(int icol = 0; icol <= irow; icol++){
+      double cov_val=0;
+      // Not in the same blocks
+      if(adapt_block_matrix_indices[icol] == block){
+        // Calculate Covariance for block
+        // https://projecteuclid.org/journals/bernoulli/volume-7/issue-2/An-adaptive-Metropolis-algorithm/bj/1080222083.full
+        cov_val = (*adaptive_covariance)(irow, icol)*Npars/5.6644;
+        cov_val += par_means_prev[irow]*par_means_prev[icol]; //First we remove the current means
+        cov_val = (cov_val*steps_post_burn+_fCurrVal[irow]*_fCurrVal[icol])/(steps_post_burn+1); //Now get mean(iRow*iCol)
+        cov_val -= par_means[icol]*par_means[irow];
+        cov_val*=5.6644/Npars;
+      }
+      (*adaptive_covariance)(icol, irow) = cov_val;
+      (*adaptive_covariance)(irow, icol) = cov_val;
+    }
+  }
+}
+
+// ********************************************
+void AdaptiveMCMCHandler::Print() {
+// ********************************************
+  MACH3LOG_INFO("Adaptive MCMC Info:");
+  MACH3LOG_INFO("Throwing from New Matrix from Step : {}", start_adaptive_throw);
+  MACH3LOG_INFO("Adaption Matrix Start Update       : {}", start_adaptive_update);
+  MACH3LOG_INFO("Adaption Matrix Ending Updates     : {}", end_adaptive_update);
+  MACH3LOG_INFO("Steps Between Updates              : {}", adaptive_update_step);
+}
+
+
+} //end adaptive_mcmc

covariance/AdaptiveMCMCHandler.h

@@ -0,0 +1,84 @@
+#pragma once
+
+// MaCh3 Includes
+#include "manager/MaCh3Logger.h"
+#include "manager/manager.h"
+#include "covariance/CovarianceUtils.h"
+
+namespace adaptive_mcmc{
+
+/// @brief Contains information about adaptive covariance matrix
+/// @see An adaptive Metropolis algorithm, H.Haario et al., 2001 for more info!
+///@details struct encapsulating all adaptive MCMC information
+class AdaptiveMCMCHandler{
+ public:
+
+  /// @brief Constructor
+  AdaptiveMCMCHandler();
+
+  /// @brief Destructor
+  ~AdaptiveMCMCHandler();
+
+  /// @brief Print all class members
+  void Print();
+
+  /// @brief Read initial values from config file
+  /// @param adapt_manager Config file from which we update matrix
+  void InitFromConfig(const YAML::Node& adapt_manager, const std::string& matrix_name_str, const int Npars);
+
+  /// @brief If we don't have a covariance matrix to start from for adaptive tune we need to make one!
+  void CreateNewAdaptiveCovariance(const int Npars);
+
+  /// @brief HW: sets adaptive block matrix
+  /// @param block_indices Values for sub-matrix blocks
+  void SetAdaptiveBlocks(std::vector<std::vector<int>> block_indices, const int Npars);
+
+  /// @brief HW: Save adaptive throw matrix to file
+  void SaveAdaptiveToFile(const TString& outFileName, const TString& systematicName);
+
+  /// @brief sets throw matrix from a file
+  /// @param matrix_file_name name of file matrix lives in
+  /// @param matrix_name name of matrix in file
+  /// @param means_name name of means vec in file
+  void SetThrowMatrixFromFile(const std::string& matrix_file_name,
+                              const std::string& matrix_name,
+                              const std::string& means_name,
+                              bool& use_adaptive,
+                              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
+  /// @param _fCurrVal Value of each parameter necessary for updating throw matrix
+  void UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int steps_post_burn, const int Npars);
+
+  /// Meta variables related to adaption run time
+  /// When do we start throwing
+  int start_adaptive_throw;
+
+  /// When do we stop update the adaptive matrix
+  int start_adaptive_update;
+
+  /// Steps between changing throw matrix
+  int end_adaptive_update;
+
+  /// Steps between changing throw matrix
+  int adaptive_update_step;
+
+  /// Indices for block-matrix adaption
+  std::vector<int> adapt_block_matrix_indices;
+
+  /// Size of blocks for adaption
+  std::vector<int> adapt_block_sizes;
+
+  // Variables directedly linked to adaption
+  /// Mean values for all parameters
+  std::vector<double> par_means;
+
+  /// Full adaptive covariance matrix
+  TMatrixDSym* adaptive_covariance;
+};
+
+} // adaptive_mcmc namespace

covariance/CMakeLists.txt

@@ -4,14 +4,15 @@ set(HEADERS
     covarianceBase.h
     covarianceOsc.h
     covarianceXsec.h
-    adaptiveMCMCStruct.h
+    AdaptiveMCMCHandler.h
 )
 
 add_library(Covariance SHARED
     ThrowParms.cpp
     covarianceBase.cpp
     covarianceOsc.cpp
     covarianceXsec.cpp
+    AdaptiveMCMCHandler.cpp
 )
 
 set_target_properties(Covariance PROPERTIES

covariance/CovarianceUtils.h

@@ -5,23 +5,19 @@
 #include "TMatrixDSym.h"
 #include "TVectorT.h"
 #include "TVectorD.h"
-#include "TCanvas.h"
 #include "TH1D.h"
 #include "TH2D.h"
 #include "TTree.h"
 #include "TFile.h"
-#include "TAxis.h"
 #include "TRandom3.h"
 #include "TMath.h"
-#include "math.h"
 #include "TDecompChol.h"
 #include "TStopwatch.h"
 #include "TMatrix.h"
 #include "TMatrixDSymEigen.h"
 #include "TMatrixDEigen.h"
 #include "TDecompSVD.h"
 
-
 #ifdef MULTITHREAD
 #include "omp.h"
 #endif