Skip to content

Commit

Permalink
Merge pull request #129 from mach3-software/feature_AdaptionHandlerUp…
Browse files Browse the repository at this point in the history 
…date

Tidy Cov handlers
  • Loading branch information
@KSkwarczynski
KSkwarczynski authored Sep 20, 2024
2 parents 1bd5a05 + 495b6ee commit 70375bd
Show file tree
Hide file tree
Showing 7 changed files with 121 additions and 83 deletions.
2 changes: 1 addition & 1 deletion cmake/Modules/MaCh3Dependencies.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ else()
CPMFindPackage(
NAME CUDAProb3
GITHUB_REPOSITORY "mach3-software/CUDAProb3"
GIT_TAG "feature_cleanup"
GIT_TAG "main"
DOWNLOAD_ONLY YES
)
LIST(APPEND MaCh3_Oscillator_ENABLED "CUDAProb3")
Expand Down
55 changes: 53 additions & 2 deletions covariance/AdaptiveMCMCHandler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ AdaptiveMCMCHandler::AdaptiveMCMCHandler() {
start_adaptive_update = 0;
end_adaptive_update = 1;
adaptive_update_step = 1000;
total_steps = 0;

par_means = {};
adaptive_covariance = nullptr;
Expand All @@ -25,7 +26,25 @@ AdaptiveMCMCHandler::~AdaptiveMCMCHandler() {
// ********************************************
bool AdaptiveMCMCHandler::InitFromConfig(const YAML::Node& adapt_manager, const std::string& matrix_name_str, const int Npars) {
// ********************************************
// setAdaptionDefaults();
/*
* HW: Idea is that adaption can simply read the YAML config
* Options :
* External Info:
* UseExternalMatrix [bool] : Use an external matrix
* ExternalMatrixFileName [str] : Name of file containing external info
* ExternalMatrixName [str] : Name of external Matrix
* ExternalMeansName [str] : Name of external means vector [for updates]
*
* General Info:
* DoAdaption [bool] : Do we want to do adaption?
* AdaptionStartThrow [int] : Step we start throwing adaptive matrix from
* AdaptionEndUpdate [int] : Step we stop updating adaptive matrix
* AdaptionStartUpdate [int] : Do we skip the first N steps?
* AdaptionUpdateStep [int] : Number of steps between matrix updates
* Adaption blocks [vector<vector<int>>] : Splits the throw matrix into several block matrices
*/

// setAdaptionDefaults();
if(!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 false;
Expand Down Expand Up @@ -176,10 +195,12 @@ void AdaptiveMCMCHandler::SetThrowMatrixFromFile(const std::string& matrix_file_
}

// ********************************************
void AdaptiveMCMCHandler::UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int steps_post_burn, const int Npars) {
void AdaptiveMCMCHandler::UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int Npars) {
// ********************************************
std::vector<double> par_means_prev = par_means;

int steps_post_burn = total_steps - start_adaptive_update;

#ifdef MULTITHREAD
#pragma omp parallel for
#endif
Expand Down Expand Up @@ -212,6 +233,36 @@ void AdaptiveMCMCHandler::UpdateAdaptiveCovariance(const std::vector<double>& _f
}
}

// ********************************************
bool AdaptiveMCMCHandler::IndivStepScaleAdapt() {
// ********************************************
if(total_steps == start_adaptive_throw) return true;
else return false;
}

// ********************************************
bool AdaptiveMCMCHandler::UpdateMatrixAdapt() {
// ********************************************
if(total_steps >= start_adaptive_throw &&
(total_steps - start_adaptive_throw)%adaptive_update_step == 0) return true;
else return false;
}

// ********************************************
bool AdaptiveMCMCHandler::SkipAdaption() {
// ********************************************
if(total_steps > end_adaptive_update ||
total_steps< start_adaptive_update) return true;
else return false;
}

// ********************************************
bool AdaptiveMCMCHandler::AdaptionUpdate() {
// ********************************************
if(total_steps <= start_adaptive_throw) return true;
else return false;
}

// ********************************************
void AdaptiveMCMCHandler::Print() {
// ********************************************
Expand Down
20 changes: 16 additions & 4 deletions covariance/AdaptiveMCMCHandler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,22 @@ class AdaptiveMCMCHandler{
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);
void UpdateAdaptiveCovariance(const std::vector<double>& _fCurrVal, const int Npars);

/// @brief Tell whether we want reset step scale or not
bool IndivStepScaleAdapt();

/// @brief Tell whether matrix should be updated
bool UpdateMatrixAdapt();

/// @brief To be fair not a clue...
bool AdaptionUpdate();

/// @brief Tell if we are Skipping Adaption
bool SkipAdaption();

/// Meta variables related to adaption run time
/// When do we start throwing
Expand All @@ -78,6 +87,9 @@ class AdaptiveMCMCHandler{

/// Full adaptive covariance matrix
TMatrixDSym* adaptive_covariance;

/// Total number of MCMC steps
int total_steps;
};

} // adaptive_mcmc namespace
20 changes: 19 additions & 1 deletion covariance/PCAHandler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,28 @@ PCAHandler::~PCAHandler() {
// ********************************************
void PCAHandler::ConstructPCA(TMatrixDSym * covMatrix, const int firstPCAd, const int lastPCAd, const double eigen_thresh, int& _fNumParPCA) {
// ********************************************

FirstPCAdpar = firstPCAd;
LastPCAdpar = lastPCAd;
eigen_threshold = eigen_thresh;

// Check that covariance matrix exists
if (covMatrix == NULL) {
MACH3LOG_ERROR("Covariance matrix for has not yet been set");
MACH3LOG_ERROR("Can not construct PCA until it is set");
throw MaCh3Exception(__FILE__ , __LINE__ );
}

if(FirstPCAdpar > covMatrix->GetNrows()-1 || LastPCAdpar>covMatrix->GetNrows()-1) {
MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are higher than the number of parameters");
MACH3LOG_ERROR("first: {} last: {}, params: {}", FirstPCAdpar, LastPCAdpar, covMatrix->GetNrows()-1);
throw MaCh3Exception(__FILE__ , __LINE__ );
}
if(FirstPCAdpar < 0 || LastPCAdpar < 0){
MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are less than 0 but not default -999");
MACH3LOG_ERROR("first: {} last: {}", FirstPCAdpar, LastPCAdpar);
throw MaCh3Exception(__FILE__ , __LINE__ );
}

MACH3LOG_INFO("PCAing parameters {} through {} inclusive", FirstPCAdpar, LastPCAdpar);
int numunpcadpars = covMatrix->GetNrows()-(LastPCAdpar-FirstPCAdpar+1);

Expand Down
2 changes: 1 addition & 1 deletion covariance/PCAHandler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#endif

/// @brief Class responsible for handling Principal Component Analysis (PCA) of covariance matrix
/// @see For more details, visit the [Wiki](https://github.com/mach3-software/MaCh3/wiki/03.-Eigen-Decomposition-%E2%80%90-PCA).
class PCAHandler{
public:

Expand All @@ -35,7 +36,6 @@ class PCAHandler{
~PCAHandler();

/// @brief CW: Calculate eigen values, prepare transition matrices and remove param based on defined threshold
/// @see For more details, visit the [Wiki](https://github.com/mach3-software/MaCh3/wiki/03.-Eigen-Decomposition-%E2%80%90-PCA).
void ConstructPCA(TMatrixDSym * covMatrix, const int firstPCAd, const int lastPCAd, const double eigen_thresh, int& _fNumParPCA);

#ifdef DEBUG_PCA
Expand Down
70 changes: 11 additions & 59 deletions covariance/covarianceBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,13 +117,6 @@ covarianceBase::~covarianceBase(){
void covarianceBase::ConstructPCA() {
// ********************************************

// Check that covariance matrix exists
if (covMatrix == NULL) {
MACH3LOG_ERROR("Covariance matrix for {} has not yet been set", matrixName);
MACH3LOG_ERROR("Can not construct PCA until it is set");
throw MaCh3Exception(__FILE__ , __LINE__ );
}

//Check whether first and last pcadpar are set and if not just PCA everything
if(FirstPCAdpar == -999 || LastPCAdpar == -999){
if(FirstPCAdpar == -999 && LastPCAdpar == -999){
Expand All @@ -135,16 +128,6 @@ void covarianceBase::ConstructPCA() {
throw MaCh3Exception(__FILE__ , __LINE__ );
}
}
if(FirstPCAdpar > covMatrix->GetNrows()-1 || LastPCAdpar>covMatrix->GetNrows()-1){
MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are higher than the number of parameters");
MACH3LOG_ERROR("first: {} last: {}, params: {}", FirstPCAdpar, LastPCAdpar, covMatrix->GetNrows()-1);
throw MaCh3Exception(__FILE__ , __LINE__ );
}
if(FirstPCAdpar < 0 || LastPCAdpar < 0){
MACH3LOG_ERROR("FirstPCAdpar and LastPCAdpar are less than 0 but not default -999");
MACH3LOG_ERROR("first: {} last: {}", FirstPCAdpar, LastPCAdpar);
throw MaCh3Exception(__FILE__ , __LINE__ );
}

PCAObj.ConstructPCA(covMatrix, FirstPCAdpar, LastPCAdpar, eigen_threshold, _fNumParPCA);
// Make a note that we have now done PCA
Expand Down Expand Up @@ -199,7 +182,7 @@ void covarianceBase::init(const char *name, const char *file) {
}

// Not using adaptive by default
setAdaptionDefaults();
use_adaptive = false;
// Set the covariance matrix
size = CovMat->GetNrows();
_fNumPar = size;
Expand Down Expand Up @@ -267,7 +250,7 @@ void covarianceBase::init(const std::vector<std::string>& YAMLFile) {
_fNumPar = _fYAMLDoc["Systematics"].size();
size = _fNumPar;

setAdaptionDefaults();
use_adaptive = false;

InvertCovMatrix = new double*[_fNumPar]();
throwMatrixCholDecomp = new double*[_fNumPar]();
Expand Down Expand Up @@ -401,7 +384,7 @@ void covarianceBase::init(TMatrixDSym* covMat) {
}
}

setAdaptionDefaults();
use_adaptive = false;
setCovMatrix(covMat);

ReserveMemory(_fNumPar);
Expand Down Expand Up @@ -649,10 +632,7 @@ void covarianceBase::proposeStep() {
// Make the random numbers for the step proposal
randomize();
CorrelateSteps();
if(use_adaptive){
updateAdaptiveCovariance();
total_steps++;
}
if(use_adaptive) updateAdaptiveCovariance();
}

// ************************************************
Expand Down Expand Up @@ -1191,7 +1171,7 @@ void covarianceBase::MakePosDef(TMatrixDSym *cov) {
MACH3LOG_ERROR("This indicates that something is wrong with the input matrix");
throw MaCh3Exception(__FILE__ , __LINE__ );
}
if(total_steps < 2) {
if(AdaptiveHandler.total_steps < 2) {
MACH3LOG_INFO("Had to shift diagonal {} time(s) to allow the covariance matrix to be decomposed", iAttempt);
}
//DB Resetting warning level
Expand Down Expand Up @@ -1228,7 +1208,7 @@ void covarianceBase::setThrowMatrix(TMatrixDSym *cov){
}

throwMatrix = (TMatrixDSym*)cov->Clone();
if(use_adaptive && total_steps <= AdaptiveHandler.start_adaptive_throw) makeClosestPosDef(throwMatrix);
if(use_adaptive && AdaptiveHandler.AdaptionUpdate()) makeClosestPosDef(throwMatrix);
else MakePosDef(throwMatrix);

TDecompChol TDecompChol_throwMatrix(*throwMatrix);
Expand Down Expand Up @@ -1264,37 +1244,10 @@ void covarianceBase::updateThrowMatrix(TMatrixDSym *cov){
setThrowMatrix(cov);
}


// ********************************************
void covarianceBase::setAdaptionDefaults(){
// ********************************************
// Puts adaptive MCMC default attributes somewhere obvious
use_adaptive = false;
total_steps = 0;
}

// ********************************************
// HW : Here be adaption
void covarianceBase::initialiseAdaption(const YAML::Node& adapt_manager){
// ********************************************
/*
HW: Idea is that adaption can simply read the YAML config
Options :
External Info:
* UseExternalMatrix [bool] : Use an external matrix
* ExternalMatrixFileName [str] : Name of file containing external info
* ExternalMatrixName [str] : Name of external Matrix
* ExternalMeansName [str] : Name of external means vector [for updates]
General Info:
* DoAdaption [bool] : Do we want to do adaption?
* AdaptionStartThrow [int] : Step we start throwing adaptive matrix from
* AdaptionEndUpdate [int] : Step we stop updating adaptive matrix
* AdaptionStartUpdate [int] : Do we skip the first N steps?
* AdaptionUpdateStep [int] : Number of steps between matrix updates
* Adaption blocks [vector<vector<int>>] : Splits the throw matrix into several block matrices
*/

// need to cast matrixName to string
std::string matrix_name_str(matrixName);

Expand Down Expand Up @@ -1331,21 +1284,21 @@ void covarianceBase::updateAdaptiveCovariance(){
// First we update the total means

// Skip this if we're at a large number of steps
if(total_steps>AdaptiveHandler.end_adaptive_update || total_steps<AdaptiveHandler.start_adaptive_update) return;
if(AdaptiveHandler.SkipAdaption()) return;

int steps_post_burn = total_steps - AdaptiveHandler.start_adaptive_update;
// Call main adaption function
AdaptiveHandler.UpdateAdaptiveCovariance(_fCurrVal, steps_post_burn, _fNumPar);
AdaptiveHandler.UpdateAdaptiveCovariance(_fCurrVal, _fNumPar);

//This is likely going to be the slow bit!
if(total_steps == AdaptiveHandler.start_adaptive_throw) {
if(AdaptiveHandler.IndivStepScaleAdapt()) {
resetIndivStepScale();
}

if(total_steps >= AdaptiveHandler.start_adaptive_throw && (total_steps-AdaptiveHandler.start_adaptive_throw)%AdaptiveHandler.adaptive_update_step==0) {
if(AdaptiveHandler.UpdateMatrixAdapt()) {
TMatrixDSym* update_matrix = static_cast<TMatrixDSym*>(AdaptiveHandler.adaptive_covariance->Clone());
updateThrowMatrix(update_matrix); //Now we update and continue!
}
AdaptiveHandler.total_steps++;
}

// ********************************************
Expand Down Expand Up @@ -1441,7 +1394,6 @@ TH2D* covarianceBase::GetCorrelationMatrix() {
// KS: After step scale, prefit etc. value were modified save this modified config.
void covarianceBase::SaveUpdatedMatrixConfig() {
// ********************************************

if (!_fYAMLDoc)
{
MACH3LOG_CRITICAL("Yaml node hasn't been initialised for matrix {}, something is not right", matrixName);
Expand Down
Loading

0 comments on commit 70375bd

Please sign in to comment.