diff --git a/covariance/CMakeLists.txt b/covariance/CMakeLists.txt index 5664f44c..60a171bb 100644 --- a/covariance/CMakeLists.txt +++ b/covariance/CMakeLists.txt @@ -4,6 +4,7 @@ set(HEADERS covarianceBase.h covarianceOsc.h covarianceXsec.h + adaptiveMCMCStruct.h ) add_library(Covariance SHARED diff --git a/covariance/adaptiveMCMCStruct.h b/covariance/adaptiveMCMCStruct.h new file mode 100644 index 00000000..02ef6d36 --- /dev/null +++ b/covariance/adaptiveMCMCStruct.h @@ -0,0 +1,54 @@ +#pragma once + +// MaCh3 Includes +#include "manager/MaCh3Logger.h" + +// ROOT Includes +#include "TMatrixDSym.h" + +// C++ Includes +#include + +/// @brief Contains information about adaptive covariance matrix +/// @see An adaptive Metropolis algorithm, H.Haario et al., 2001 for more info! + + +namespace adaptive_mcmc{ + +///@brief struct encapsulating all adaptive MCMC information +struct AdaptiveMCMCStruct{ + // 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 adapt_block_matrix_indices; + ///Size of blocks for adaption + std::vector adapt_block_sizes; + + // Variables directely linked to adaption + /// Mean values for all parameters + std::vector par_means; + /// Full adaptive covariance matrix + TMatrixDSym* adaptive_covariance; +}; + + +///@brief Basic printer method for AMCMC struct +inline void print_adaptive_struct(AdaptiveMCMCStruct adaptive_struct){ + MACH3LOG_INFO("Adaptive MCMC Info:"); + MACH3LOG_INFO("Throwing from New Matrix from Step : {}", adaptive_struct.start_adaptive_throw); + MACH3LOG_INFO("Adaption Matrix Start Update : {}", adaptive_struct.start_adaptive_update); + MACH3LOG_INFO("Adaption Matrix Ending Updates : {}", adaptive_struct.end_adaptive_update); + MACH3LOG_INFO("Steps Between Updates : {}", adaptive_struct.adaptive_update_step); +} + + + + +} // adaptive_mcmc namespace diff --git a/covariance/covarianceBase.cpp b/covariance/covarianceBase.cpp index 7ccf7941..60f5820f 100644 --- a/covariance/covarianceBase.cpp +++ b/covariance/covarianceBase.cpp @@ -1,5 +1,4 @@ #include "covarianceBase.h" - // ******************************************** covarianceBase::covarianceBase(const char *name, const char *file) : inputFile(std::string(file)), pca(false) { // ******************************************** @@ -112,6 +111,11 @@ covarianceBase::~covarianceBase(){ for (int iThread = 0;iThread < nThreads; iThread++) delete random_number[iThread]; delete[] random_number; if (throwMatrix != NULL) delete throwMatrix; + + if(adaption_struct.adaptive_covariance != nullptr){ + delete adaption_struct.adaptive_covariance; + } + } // ******************************************** @@ -280,11 +284,7 @@ void covarianceBase::init(const char *name, const char *file) { } // Not using adaptive by default - use_adaptive = false; - total_steps = 0; - lower_adapt = -999; - upper_adapt = -999; - + setAdaptionDefaults(); // Set the covariance matrix size = CovMat->GetNrows(); _fNumPar = size; @@ -349,16 +349,13 @@ void covarianceBase::init(const std::vector& YAMLFile) { PrintLength = 35; - // Not using adaptive by default - use_adaptive = false; - total_steps = 0; - lower_adapt = -999; - upper_adapt = -999; // Set the covariance matrix _fNumPar = _fYAMLDoc["Systematics"].size(); size = _fNumPar; - + + setAdaptionDefaults(); + InvertCovMatrix = new double*[_fNumPar](); throwMatrixCholDecomp = new double*[_fNumPar](); // Set the defaults to true @@ -491,7 +488,8 @@ void covarianceBase::init(TMatrixDSym* covMat) { throwMatrixCholDecomp[i][j] = 0.; } } - + + setAdaptionDefaults(); setCovMatrix(covMat); ReserveMemory(_fNumPar); @@ -787,7 +785,10 @@ void covarianceBase::proposeStep() { // Make the random numbers for the step proposal randomize(); CorrelateSteps(); - if(use_adaptive && total_steps < upper_adapt) updateAdaptiveCovariance(); + if(use_adaptive){ + updateAdaptiveCovariance(); + total_steps++; + } } // ************************************************ @@ -1358,7 +1359,7 @@ void covarianceBase::setThrowMatrix(TMatrixDSym *cov){ } throwMatrix = (TMatrixDSym*)cov->Clone(); - if(total_steps <= lower_adapt) makeClosestPosDef(throwMatrix); + if(use_adaptive && total_steps <= adaption_struct.start_adaptive_throw) makeClosestPosDef(throwMatrix); else MakePosDef(throwMatrix); TDecompChol TDecompChol_throwMatrix(*throwMatrix); @@ -1393,125 +1394,276 @@ void covarianceBase::updateThrowMatrix(TMatrixDSym *cov){ setThrowMatrix(cov); } -// The setter -void covarianceBase::useSeparateThrowMatrix(TString throwMatrixFileName, TString throwMatrixName, TString meansVectorName){ -// Firstly let's check if the file exists - TFile* throwMatrixFile = new TFile(throwMatrixFileName); - resetIndivStepScale(); - if(throwMatrixFile->IsZombie()) { - MACH3LOG_ERROR("Couldn't find throw Matrix file : {}", throwMatrixFileName); - std::cerr<<__FILE__<<" : "<<__LINE__<IsZombie()){ + MACH3LOG_ERROR("Couldn't find {}", outFileName); throw; - } //We're done for now + } + TVectorD* outMeanVec = new TVectorD((int)adaption_struct.par_means.size()); + for(int i = 0; i < (int)adaption_struct.par_means.size(); i++){ + (*outMeanVec)(i)=adaption_struct.par_means[i]; + } + outFile->cd(); + adaption_struct.adaptive_covariance->Write(systematicName+"_postfit_matrix"); + outMeanVec->Write(systematicName+"_mean_vec"); + outFile->Close(); + delete outMeanVec; + delete outFile; +} - TMatrixDSym* tmp_throwMatrix = (TMatrixDSym*)throwMatrixFile->Get(throwMatrixName); - TVectorD* tmp_meansvec = (TVectorD*)throwMatrixFile->Get(meansVectorName); - if(!tmp_throwMatrix){ - std::cerr<<"ERROR : Couldn't find throw matrix "<>] : Splits the throw matrix into several block matrices + */ + + // need to cast matrixName to string + std::string matrix_name_str(matrixName); + + // setAdaptionDefaults(); + + if(GetFromManager(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(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["DoAdaption"], false)) { + MACH3LOG_INFO("Not using adaption for {}", matrix_name_str); + return; } - if(!tmp_meansvec){ - std::cerr<<"ERROR : Couldn't find means vector "<(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionStartThrow"], 10); + adaption_struct.start_adaptive_update = GetFromManager(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionStartUpdate"], 0); + adaption_struct.end_adaptive_update = GetFromManager(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionEndUpdate"], 10000); + adaption_struct.adaptive_update_step = GetFromManager(adapt_manager["AdaptionOptions"]["Settings"]["AdaptionUpdateStep"], 100); + + + adaptive_mcmc::print_adaptive_struct(adaption_struct); + + // 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>>(adapt_manager["AdaptionOptions"]["Settings"][matrix_name_str]["AdaptionUpdateStep"], {{}}); + + setAdaptiveBlocks(matrix_blocks); + + // Next let"s check for external matrices + // We"re going to grab this info from the YAML manager + if(!GetFromManager(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["UseExternalMatrix"], false)) { + MACH3LOG_INFO("Not using external matrix for {}, initialising adaption from scratch", matrix_name_str); + createNewAdaptiveCovariance(); + return; + } + + // Finally, we accept that we want to read the matrix from a file! + // + std::string external_file_name = GetFromManager(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMatrixFileName"], ""); + std::string external_matrix_name = GetFromManager(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMatrixName"], ""); + std::string external_mean_name = GetFromManager(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["ExternalMeansName"], ""); + + setThrowMatrixFromFile(external_file_name, external_matrix_name, external_mean_name); + +} + +// HW : I would like this to be less painful to use! +// First things first we need setters +void covarianceBase::setThrowMatrixFromFile(std::string matrix_file_name, std::string matrix_name, std::string means_name){ + // Lets you set the throw matrix externally + // Open file + std::unique_ptrmatrix_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 + adaption_struct.adaptive_covariance = static_cast(matrix_file->Get(matrix_name.c_str())); + if(!adaption_struct.adaptive_covariance){ + MACH3LOG_ERROR("Coukdn't find {} in {}", matrix_name, matrix_file_name); throw; } - adaptiveCovariance = (TMatrixDSym*)tmp_throwMatrix->Clone(); - par_means.resize(tmp_meansvec->GetNrows()); - for(int iMean = 0; iMean < tmp_meansvec->GetNrows(); iMean++){ - par_means[iMean] = (*tmp_meansvec)(iMean); - updateThrowMatrix(adaptiveCovariance); + MACH3LOG_INFO("Succesfully Set External Throw Matrix Stored in {}", matrix_file_name); + + setThrowMatrix(adaption_struct.adaptive_covariance); + + // Finally we grab the means vector + TVectorD* means_vector = static_cast(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(), size); + throw MaCh3Exception(__FILE__, __LINE__); + } + + adaption_struct.par_means = std::vector(size); + for(int i=0; iClose(); - delete throwMatrixFile; // Just in case + // 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(); + + std::cout<<"Set up matrix from external file"<> block_indices){ + /* + 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 + adaption_struct.adapt_block_matrix_indices = std::vector(getNpars(), 0); + + // Should also make a matrix of block sizes + adaption_struct.adapt_block_sizes = std::vector((int)block_indices.size()+1, 0); + adaption_struct.adapt_block_sizes[0]=getNpars(); + + 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<getNpars() || block_ub>getNpars()){ + MACH3LOG_ERROR("Cannot set matrix block with edges {}, {} for matrix of size {}", + block_lb, block_ub, getNpars()); + throw MaCh3Exception(__FILE__, __LINE__);; + } + + for(int ipar=block_lb; iparZero(); + adaption_struct.par_means = std::vector(size, 0); } -//HW: Truly adaptive MCMC! + + +// Truely adaptive MCMC! void covarianceBase::updateAdaptiveCovariance(){ // https://projecteuclid.org/journals/bernoulli/volume-7/issue-2/An-adaptive-Metropolis-algorithm/bj/1080222083.full // Updates adaptive matrix // First we update the total means -#ifdef MULTITHREAD -#pragma omp parallel for -#endif - for(int iRow = 0; iRow < _fNumPar; ++iRow){ - par_means_prev[iRow] = par_means[iRow]; - par_means[iRow] = (_fCurrVal[iRow]+par_means[iRow]*total_steps)/(total_steps+1); + // Skip this if we're at a large number of steps + if(total_steps>adaption_struct.end_adaptive_update || total_steps par_means_prev = adaption_struct.par_means; + + #ifdef MULTITHREAD + #pragma omp parallel for + #endif + for(int iRow=0; iRow= lower_adapt) { - updateThrowMatrix(adaptiveCovariance); //Now we update and continue! + if(total_steps>=adaption_struct.start_adaptive_throw && (total_steps-adaption_struct.start_adaptive_throw)%adaption_struct.adaptive_update_step==0) { + TMatrixDSym* update_matrix = static_cast(adaption_struct.adaptive_covariance->Clone()); + updateThrowMatrix(update_matrix); //Now we update and continue! } - // if(total_steps%1000==0 && total_steps>1000){ - // suboptimality_vals.push_back(calculatesuboptimality(covSqrt, adaptiveCovariance)); - // } -} - -void covarianceBase::initialiseNewAdaptiveChain(){ - // If we don't have a covariance matrix to start from for adaptive tune we need to make one! - adaptiveCovariance = new TMatrixDSym(_fNumPar); - par_means.resize(_fNumPar); - par_means_prev.resize(_fNumPar); - for(int i = 0; i < _fNumPar; ++i){ - par_means[i] = 0.; - par_means_prev[i] = 0.; - for(int j = 0; j <= i; ++j){ - (*adaptiveCovariance)(i,j) = 0.0; // Just make an empty matrix - (*adaptiveCovariance)(j,i) = 0.0; // Just make an empty matrix - } - } -} - -void covarianceBase::saveAdaptiveToFile(TString outFileName, 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(); - adaptiveCovariance->Write(systematicName+"_postfit_matrix"); - outMeanVec->Write(systematicName+"_mean_vec"); - outFile->Close(); - delete outFile; } //HW: Finds closest possible positive definite matrix in Frobenius Norm ||.||_frob diff --git a/covariance/covarianceBase.h b/covariance/covarianceBase.h index aa131137..6962bbc2 100644 --- a/covariance/covarianceBase.h +++ b/covariance/covarianceBase.h @@ -5,6 +5,7 @@ #include "covariance/CovarianceUtils.h" #include "covariance/ThrowParms.h" #include "manager/manager.h" +#include "covariance/adaptiveMCMCStruct.h" #ifndef _LARGE_LOGL_ /// Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should eb removed @@ -162,22 +163,21 @@ class covarianceBase { // Adaptive Step Tuning Stuff void resetIndivStepScale(); - void useSeparateThrowMatrix(TString throwMatrixName, TString throwMatrixFile, TString parameterMeansName=""); - void useSeparateThrowMatrix(); + + void initialiseAdaption(YAML::Node& adapt_manager); void saveAdaptiveToFile(TString outFileName, TString systematicName); + bool getDoAdaption(){return use_adaptive;} void setThrowMatrix(TMatrixDSym *cov); void updateThrowMatrix(TMatrixDSym *cov); inline void setNumberOfSteps(const int nsteps) { total_steps = nsteps; - if(total_steps >= lower_adapt) resetIndivStepScale(); + if(total_steps >= adaption_struct.start_adaptive_throw) resetIndivStepScale(); } - /// @brief HW: Set thresholds for MCMC steps - inline void setAdaptiveThresholds(const int low_threshold = 10000, const int up_threshold = 1000000){lower_adapt=low_threshold; upper_adapt = up_threshold;} inline TMatrixDSym *getThrowMatrix(){return throwMatrix;} inline TMatrixD *getThrowMatrix_CholDecomp(){return throwMatrix_CholDecomp;} - inline std::vector getParameterMeans(){return par_means;} + inline std::vector getParameterMeans(){return adaption_struct.par_means;} /// @brief KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting TH2D* GetCorrelationMatrix(); @@ -369,18 +369,8 @@ class covarianceBase { /// @brief KS: Custom function to perform multiplication of matrix and single element which is thread safe inline double MatrixVectorMultiSingle(double** _restrict_ matrix, const double* _restrict_ vector, const int Length, const int i); - /// @brief HW: Turn on/off true adaptive MCMC, Also set thresholds for use (having a lower threshold gives us some data to adapt from!) - void enableAdaptiveMCMC(bool enable = true){ - use_adaptive = enable; - total_steps = 0; //Set these to default values - lower_adapt = 10000; - upper_adapt = 10000000; - } - /// @brief HW: Update throw matrix used for proposal - void updateAdaptiveCovariance(); +protected: - protected: - /// @brief Initialisation of the class using matrix stored in the ROOT file void init(const char *name, const char *file); /// @brief Initialisation of the class using config /// @param YAMLFile A vector of strings representing the YAML files used for initialisation of matrix @@ -506,16 +496,28 @@ class covarianceBase { /// Throw matrix that is being used in the fit, much faster as TMatrixDSym cache miss double **throwMatrixCholDecomp; - /// @brief HW: Truly Adaptive Stuff - void initialiseNewAdaptiveChain(); - // TMatrixD* getMatrixSqrt(TMatrixDSym* inputMatrix); - // double calculateSubmodality(TMatrixD* sqrtVectorCov, TMatrixDSym* throwCov); - /// HW: Do we use adaptive MCMC or not + /// @brief sets default values for adaptive MCMC parameters + void setAdaptionDefaults(); + + /// @brief sets adaptive block matrix + /// @param block_indices Values for sub-matrix blocks + void setAdaptiveBlocks(std::vector> block_indices); + + /// @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(std::string matrix_file_name, std::string matrix_name, std::string means_name); + + /// @brief Method to update adaptive MCMC + void updateAdaptiveCovariance(); + /// @brief method to create new throw matrix + void createNewAdaptiveCovariance(); + + /// Are we using AMCMC? bool use_adaptive; int total_steps; - int lower_adapt, upper_adapt; //Thresholds for when to turn on/off adaptive MCMC - // TMatrixD* covSqrt; - std::vector par_means; - std::vector par_means_prev; - TMatrixDSym* adaptiveCovariance; + + /// Struct containing information about adaption + adaptive_mcmc::AdaptiveMCMCStruct adaption_struct; };