Merge pull request #184 from mach3-software/feature_UpdateStartFromPr…

…evious More Flexible Start FromPrevious
mach3-software · Oct 28, 2024 · d22fd56 · d22fd56
2 parents 82add79 + 1790252
commit d22fd56
Show file tree

Hide file tree

Showing 4 changed files with 81 additions and 57 deletions.
diff --git 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<double> 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
@@ -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
@@ -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<double> 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
@@ -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";};