diff --git a/mcmc/FitterBase.cpp b/mcmc/FitterBase.cpp index 96751fa37..b5f105148 100644 --- a/mcmc/FitterBase.cpp +++ b/mcmc/FitterBase.cpp @@ -300,6 +300,76 @@ void FitterBase::addSystObj(covarianceBase * const cov) { return; } + +// ******************* +void FitterBase::StartFromPreviousFit(const std::string& FitName) { +// ******************* + MACH3LOG_INFO("Getting starting position from {}", FitName); + + TFile *infile = new TFile(FitName.c_str(), "READ"); + TTree *posts = (TTree*)infile->Get("posteriors"); + int step_val = 0; + double log_val = _LARGE_LOGL_; + posts->SetBranchAddress("step",&step_val); + posts->SetBranchAddress("LogL",&log_val); + + for (size_t s = 0; s < systematics.size(); ++s) + { + TDirectory* CovarianceFolder = (TDirectory*)infile->Get("CovarianceFolder"); + + std::string ConfigName = "Config_" + std::string(systematics[s]->getName()); + TMacro *ConfigCov = (TMacro*)(CovarianceFolder->Get(ConfigName.c_str())); + // KS: Not every covariance uses yaml, if it uses yaml make sure they are identical + if (ConfigCov != nullptr) { + // Config which was in MCMC from which we are starting + YAML::Node CovSettings = TMacroToYAML(*ConfigCov); + // Config from currently used cov object + YAML::Node ConfigCurrent = systematics[s]->GetConfig(); + + if (!compareYAMLNodes(CovSettings, ConfigCurrent)) + { + MACH3LOG_ERROR("Yaml configs in previous chain and current one are different", FitName); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + } + if(ConfigCov == nullptr) delete ConfigCov; + + std::vector branch_vals(systematics[s]->GetNumParams()); + for (int i = 0; i < systematics[s]->GetNumParams(); ++i) { + branch_vals[i] = _BAD_DOUBLE_; + posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), &branch_vals[i]); + } + posts->GetEntry(posts->GetEntries()-1); + + for (int i = 0; i < systematics[s]->GetNumParams(); ++i) { + if(branch_vals[i] == _BAD_DOUBLE_) + { + MACH3LOG_ERROR("Parameter {} is unvitalised with value {}", i, branch_vals[i]); + MACH3LOG_ERROR("Please check more precisely chain you passed {}", FitName); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + } + systematics[s]->setParameters(branch_vals); + systematics[s]->acceptStep(); + + MACH3LOG_INFO("Printing new starting values for: {}", systematics[s]->getName()); + systematics[s]->printNominalCurrProp(); + + CovarianceFolder->Close(); + delete CovarianceFolder; + } + logLCurr = log_val; + + infile->Close(); + delete infile; + + for (size_t s = 0; s < systematics.size(); ++s) { + if(systematics[s]->getDoAdaption()){ //Use separate throw matrix for xsec + systematics[s]->setNumberOfSteps(step_val); + } + } +} + // ******************* // Process the MCMC output to get postfit etc void FitterBase::ProcessMCMC() { diff --git a/mcmc/FitterBase.h b/mcmc/FitterBase.h index 14d9d365e..71d33aec0 100644 --- a/mcmc/FitterBase.h +++ b/mcmc/FitterBase.h @@ -57,6 +57,11 @@ class FitterBase { /// @warning Code uses TH2Poly void RunSigmaVar(); + /// @brief Allow to start from previous fit/chain + /// @param FitName Name of previous chain + /// @todo implement some check that number of params matches etc + virtual void StartFromPreviousFit(const std::string& FitName); + /// @brief Get name of class virtual inline std::string GetName()const {return "FitterBase";}; protected: diff --git a/mcmc/mcmc.cpp b/mcmc/mcmc.cpp index bc3c582bc..5269dc0cf 100644 --- a/mcmc/mcmc.cpp +++ b/mcmc/mcmc.cpp @@ -217,69 +217,18 @@ void mcmc::PrintProgress() { // ******************* void mcmc::StartFromPreviousFit(const std::string& FitName) { // ******************* - MACH3LOG_INFO("Getting starting position from {}", FitName); + // Use base class + FitterBase::StartFromPreviousFit(FitName); + // For MCMC we also need to set stepStart TFile *infile = new TFile(FitName.c_str(), "READ"); TTree *posts = (TTree*)infile->Get("posteriors"); int step_val = 0; - double log_val = _LARGE_LOGL_; - posts->SetBranchAddress("step",&step_val); - posts->SetBranchAddress("LogL",&log_val); - - for (size_t s = 0; s < systematics.size(); ++s) - { - TDirectory* CovarianceFolder = (TDirectory*)infile->Get("CovarianceFolder"); - - std::string ConfigName = "Config_" + std::string(systematics[s]->getName()); - TMacro *ConfigCov = (TMacro*)(CovarianceFolder->Get(ConfigName.c_str())); - // KS: Not every covariance uses yaml, if it uses yaml make sure they are identical - if (ConfigCov != nullptr) { - // Config which was in MCMC from which we are starting - YAML::Node CovSettings = TMacroToYAML(*ConfigCov); - // Config from currently used cov object - YAML::Node ConfigCurrent = systematics[s]->GetConfig(); - - if (!compareYAMLNodes(CovSettings, ConfigCurrent)) - { - MACH3LOG_ERROR("Yaml configs in previous chain and current one are different", FitName); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - } - if(ConfigCov == nullptr) delete ConfigCov; - - std::vector branch_vals(systematics[s]->GetNumParams()); - for (int i = 0; i < systematics[s]->GetNumParams(); ++i) { - branch_vals[i] = _BAD_DOUBLE_; - posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), &branch_vals[i]); - } - posts->GetEntry(posts->GetEntries()-1); - - for (int i = 0; i < systematics[s]->GetNumParams(); ++i) { - if(branch_vals[i] == _BAD_DOUBLE_) - { - MACH3LOG_ERROR("Parameter {} is unvitalised with value {}", i, branch_vals[i]); - MACH3LOG_ERROR("Please check more precisely chain you passed {}", FitName); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - } - systematics[s]->setParameters(branch_vals); - systematics[s]->acceptStep(); - MACH3LOG_INFO("Printing new starting values for: {}", systematics[s]->getName()); - systematics[s]->printNominalCurrProp(); + posts->SetBranchAddress("step",&step_val); + posts->GetEntry(posts->GetEntries()-1); - CovarianceFolder->Close(); - delete CovarianceFolder; - } stepStart = step_val; - logLCurr = log_val; - infile->Close(); delete infile; - - for (size_t s = 0; s < systematics.size(); ++s) { - if(systematics[s]->getDoAdaption()){ //Use separate throw matrix for xsec - systematics[s]->setNumberOfSteps(step_val); - } - } } diff --git a/mcmc/mcmc.h b/mcmc/mcmc.h index 223ac89e4..9bb963aca 100644 --- a/mcmc/mcmc.h +++ b/mcmc/mcmc.h @@ -23,7 +23,7 @@ class mcmc : public FitterBase { /// @brief Allow to start from previous fit/chain /// @param FitName Name of previous chain /// @todo implement some check that number of params matches etc - void StartFromPreviousFit(const std::string& FitName); + void StartFromPreviousFit(const std::string& FitName) override; /// @brief Get name of class inline std::string GetName()const {return "MCMC";};