From 7fb5b6b42f8d9fd654f5a53728b548663e80ceed Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Sat, 9 Nov 2024 23:15:52 +0000 Subject: [PATCH 1/4] Implement simple MCMC Averaging --- .mailmap | 2 + Diagnostics/DiagMCMC.cpp | 5 +-- Diagnostics/ProcessMCMC.cpp | 51 ++++++++------------- Doc/bibliography.bib | 11 +++++ mcmc/MCMCProcessor.cpp | 3 +- mcmc/MCMCProcessor.h | 9 +++- mcmc/StatisticalUtils.h | 89 +++++++++++++++++++++++++++++++++++++ 7 files changed, 130 insertions(+), 40 deletions(-) diff --git a/.mailmap b/.mailmap index 516369c14..3fecaeb98 100644 --- a/.mailmap +++ b/.mailmap @@ -19,6 +19,8 @@ Henry Wallace Henry Wallace <67589487+henry-wallace-ph Henry Wallace Henry Wallace Henry Wallace Henry Wallace Henry Wallace Henry Wallace +Henry Wallace Henry Wallace +Henry Wallace Henry Wallace # Ewan Miller Ewan Miller Ewan diff --git a/Diagnostics/DiagMCMC.cpp b/Diagnostics/DiagMCMC.cpp index f1245b733..854b36536 100644 --- a/Diagnostics/DiagMCMC.cpp +++ b/Diagnostics/DiagMCMC.cpp @@ -12,8 +12,7 @@ void DiagMCMC(const std::string& inputFile, const std::string& config) YAML::Node Settings = YAML::LoadFile(config); // Make the processor - MCMCProcessor* Processor = new MCMCProcessor(inputFile); - + auto Processor = std::make_unique(inputFile); Processor->SetOutputSuffix("_MCMC_Diag"); //KS:Turn off plotting detector and some other setting Processor->SetExcludedTypes(GetFromManager>(Settings["DiagMCMC"]["ExcludedTypes"], {""})); @@ -27,8 +26,6 @@ void DiagMCMC(const std::string& inputFile, const std::string& config) //KS: finally call main method Processor->DiagMCMC(); - - delete Processor; } int main(int argc, char *argv[]) { diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 72bb3ac1b..37b3b34df 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -16,7 +16,7 @@ inline void ReweightPrior(MCMCProcessor* Processor); /// @brief KS: Convert TMatrix to TH2D, mostly useful for making fancy plots inline TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title); /// @brief KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution -inline void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TString canvasname); +inline void KolmogorovSmirnovTest(const std::vector>& Processor, TCanvas* Posterior, TString canvasname); int nFiles; std::vector FileNames; @@ -96,7 +96,9 @@ void ProcessMCMC(const std::string& inputFile) if(Settings["Thinning"]) { bool DoThin = Settings["Thinning"][0].as(); - if(DoThin) Processor->ThinMCMC(Settings["Thinning"][1].as()); + if(DoThin){ + Processor->ThinMCMC(Settings["Thinning"][1].as(), GetFromManager(Settings["Average"], false)); + } } // Make the postfit Processor->MakePostfit(); @@ -152,16 +154,15 @@ void MultipleProcessMCMC() YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; - const Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2}; - //const Style_t PosteriorStyle[] = {kSolid, kDashed, kDotted}; + constexpr Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2}; + //constexpr Style_t PosteriorStyle[] = {kSolid, kDashed, kDotted}; nFiles = int(FileNames.size()); - MCMCProcessor** Processor; - Processor = new MCMCProcessor*[nFiles]; + std::vector> Processor(nFiles); for (int ik = 0; ik < nFiles; ik++) { MACH3LOG_INFO("File for study: {}", FileNames[ik]); // Make the processor - Processor[ik] = new MCMCProcessor(FileNames[ik]); + Processor[ik] = std::make_unique(FileNames[ik]); Processor[ik]->SetOutputSuffix(("_" + std::to_string(ik)).c_str()); Processor[ik]->SetExcludedTypes(GetFromManager>(Settings["ExcludedTypes"], {""})); @@ -212,8 +213,8 @@ void MultipleProcessMCMC() for(int i = 0; i < Processor[0]->GetNParams(); ++i) { // This holds the posterior density - TH1D **hpost = new TH1D*[nFiles]; - TLine **hpd = new TLine*[nFiles]; + std::vector hpost(nFiles); + std::vector hpd(nFiles); hpost[0] = static_cast(Processor[0]->GetHpost(i)->Clone()); bool Skip = false; @@ -235,8 +236,6 @@ void MultipleProcessMCMC() for (int ik = 0; ik < nFiles; ik++) delete hpost[ik]; - delete[] hpost; - delete[] hpd; continue; } for (int ik = 0; ik < nFiles; ik++) @@ -258,20 +257,20 @@ void MultipleProcessMCMC() Processor[0]->GetNthParameter(i, Prior, PriorError, Title); // Now make the TLine for the Asimov - TLine *Asimov = new TLine(Prior, hpost[0]->GetMinimum(), Prior, hpost[0]->GetMaximum()); + auto Asimov = std::make_unique(Prior, hpost[0]->GetMinimum(), Prior, hpost[0]->GetMaximum()); Asimov->SetLineColor(kRed-3); Asimov->SetLineWidth(2); Asimov->SetLineStyle(kDashed); // Make a nice little TLegend - TLegend *leg = new TLegend(0.12, 0.7, 0.6, 0.97); + auto leg = std::make_unique(0.12, 0.7, 0.6, 0.97); leg->SetTextSize(0.03f); leg->SetFillColor(0); leg->SetFillStyle(0); leg->SetLineColor(0); leg->SetLineStyle(0); TString asimovLeg = Form("#splitline{Prior}{x = %.2f , #sigma = %.2f}", Prior, PriorError); - leg->AddEntry(Asimov, asimovLeg, "l"); + leg->AddEntry(Asimov.get(), asimovLeg, "l"); for (int ik = 0; ik < nFiles; ik++) { @@ -296,16 +295,11 @@ void MultipleProcessMCMC() leg->Draw("same"); Posterior->cd(); Posterior->Print(canvasname); - - delete Asimov; - delete leg; for (int ik = 0; ik < nFiles; ik++) { delete hpost[ik]; delete hpd[ik]; } - delete[] hpost; - delete[] hpd; }//End loop over parameters // Finally draw the parameter plot onto the PDF @@ -321,8 +315,6 @@ void MultipleProcessMCMC() Posterior->Print(canvasname); delete Posterior; - for (int ik = 0; ik < nFiles; ik++) delete Processor[ik]; - delete[] Processor; } // KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials @@ -627,7 +619,7 @@ TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title) } //KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution -void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TString canvasname) +void KolmogorovSmirnovTest(const std::vector>& Processor, TCanvas* Posterior, TString canvasname) { const Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2}; const Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted}; @@ -702,12 +694,9 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin CumulativeDistribution[ik]->SetBinContent(NumberOfBins+1, 1.); } - int* TestStatBin = new int[nFiles]; - double* TestStatD = new double[nFiles]; - TLine **LineD = new TLine*[nFiles]; - - for (int ik = 0 ; ik < nFiles; ik++) { TestStatBin[ik] = 0; TestStatD[ik] = -999;} - + std::vector TestStatBin(nFiles, 0); + std::vector TestStatD(nFiles, -999); + std::vector LineD(nFiles); //Find KS statistic for (int ik = 1 ; ik < nFiles; ik++) { @@ -736,7 +725,7 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin for (int ik = 0 ; ik < nFiles; ik++) CumulativeDistribution[ik]->Draw("SAME"); - TLegend *leg = new TLegend(0.15, 0.7, 0.5, 0.90); + auto leg = std::make_unique(0.15, 0.7, 0.5, 0.90); leg->SetTextSize(0.04f); for (int ik = 0; ik < nFiles; ik++) leg->AddEntry(CumulativeDistribution[ik], TitleNames[ik].c_str(), "l"); @@ -755,7 +744,6 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin Posterior->cd(); Posterior->Print(canvasname); - delete leg; for (int ik = 0; ik < nFiles; ik++) { delete hpost[ik]; @@ -764,8 +752,5 @@ void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TStrin } delete[] hpost; delete[] CumulativeDistribution; - delete[] LineD; - delete[] TestStatBin; - delete[] TestStatD; } //End loop over parameter } diff --git a/Doc/bibliography.bib b/Doc/bibliography.bib index 79fac0844..35ec3846e 100755 --- a/Doc/bibliography.bib +++ b/Doc/bibliography.bib @@ -180,3 +180,14 @@ @article{2011ThinningMCMC doi = {10.1111/j.2041-210X.2011.00131.x}, url = {https://doi.org/10.1111/j.2041-210X.2011.00131.x}, } + +@article{nesterov2009primal, + author = {Nesterov, Yurii}, + title = {Primal-dual subgradient methods for convex problems}, + journal = {Mathematical Programming}, + volume = {120}, + pages = {221--259}, + year = {2009}, + doi = {10.1007/s10107-007-0149-x}, + note = {Received 29 September 2005; Accepted 13 January 2007; Published 19 June 2007; Issue Date August 2009} +} diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index 94b8e3846..46f14a73a 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -263,8 +263,7 @@ void MCMCProcessor::MakePostfit() { // nDraw is number of draws we want to do for (int i = 0; i < nDraw; ++i) { - if (i % (nDraw/5) == 0) - { + if (i % (nDraw/5) == 0) { MaCh3Utils::PrintProgressBar(i, nDraw); } OutputFile->cd(); diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 2a5a106c3..6b6143357 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -161,7 +161,14 @@ class MCMCProcessor { /// @brief Thin MCMC Chain, to save space and maintain low autocorrelations. /// @param ThinningCut every which entry you want to thin - inline void ThinMCMC(const int ThinningCut) {ThinningMCMC(MCMCFile+".root", ThinningCut);}; + /// @param Average If true will perform MCMC averaging instead of thinning + inline void ThinMCMC(const int ThinningCut, const bool Average = false) { + if(Average){ + AveragingMCMC(MCMCFile+".root", ThinningCut); + } else { + ThinningMCMC(MCMCFile+".root", ThinningCut); + } + }; /// @brief KS: Perform MCMC diagnostic including Autocorrelation, Trace etc. void DiagMCMC(); diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h index 4312b6b3a..94a9ae225 100644 --- a/mcmc/StatisticalUtils.h +++ b/mcmc/StatisticalUtils.h @@ -584,3 +584,92 @@ inline void ThinningMCMC(const std::string& FilePath, const int ThinningCut) { MACH3LOG_INFO("Thinned TTree saved and overwrote original in: {}", TempFilePath); } + + +// ******************** +/// @brief Average MCMC Chain entries, to reduce autocorrelation and improve estimate stability. +/// +/// @param FilePath Path to the MCMC chain you want to average +/// @param AveragingCut Number of consecutive entries to average +/// @cite nesterov2009primal +/// @warning Averaging is done over entries, so this may not work well for merged chains or non-sequential entries. +/// @note The resulting file overwrites the original file. +inline void AveragingMCMC(const std::string& FilePath, const int AveragingCut) { +// ******************** + if (AveragingCut <= 1) { + MACH3LOG_WARN("Invalid AveragingCut value ({}). It must be greater than 1.", AveragingCut); + throw MaCh3Exception(__FILE__, __LINE__); + } + + // Define the path for the temporary averaged file + std::string TempFilePath = "Averaged_" + FilePath; + int ret = system(("cp " + FilePath + " " + TempFilePath).c_str()); + if (ret != 0) { + MACH3LOG_WARN("Error: system call to copy file failed with code {}", ret); + } + + TFile *inFile = TFile::Open(TempFilePath.c_str(), "UPDATE"); + if (!inFile || inFile->IsZombie()) { + MACH3LOG_ERROR("Error opening file: {}", TempFilePath); + throw MaCh3Exception(__FILE__, __LINE__); + } + + TTree *inTree = inFile->Get("posteriors"); + if (!inTree) { + MACH3LOG_ERROR("Error: TTree 'posteriors' not found in file."); + inFile->ls(); + inFile->Close(); + throw MaCh3Exception(__FILE__, __LINE__); + } + + // Clone the structure without data + TTree *outTree = inTree->CloneTree(0); + Long64_t nEntries = inTree->GetEntries(); + int nBranches = inTree->GetNbranches(); + + std::vector AvgValues(nBranches, 0.0); + std::vector BranchData(nBranches, 0.0); + + // Loop over branches and set up data pointers + TObjArray *branches = inTree->GetListOfBranches(); + for (int b = 0; b < branches->GetEntries(); ++b) { + TBranch *branch = static_cast(branches->At(b)); + inTree->SetBranchAddress(branch->GetName(), &BranchData[b]); + } + + MACH3LOG_INFO("Averaging every {} entries. Total entries: {}", AveragingCut, nEntries); + + int batchCount = 0; + for (Long64_t i = 0; i < nEntries; i += AveragingCut) { + // Reset to 0 + std::fill(AvgValues.begin(), AvgValues.end(), 0.0); + + /// @todo Here we perform normal mean. But in literature there are available weighted means and more complex operations + int entriesInBatch = 0; + for (int j = 0; j < AveragingCut && (i + j) < nEntries; ++j) { + inTree->GetEntry(i + j); + for (int b = 0; b < branches->GetEntries(); ++b) { + AvgValues[b] += BranchData[b]; + } + ++entriesInBatch; + } + for (int b = 0; b < branches->GetEntries(); ++b) { + AvgValues[b] /= entriesInBatch; + BranchData[b] = AvgValues[b]; + } + + outTree->Fill(); + + if (i % (nEntries/10) == 0) { + MaCh3Utils::PrintProgressBar(batchCount, nEntries / AveragingCut); + } + ++batchCount; + } + + // Write the averaged tree back to the file + inFile->WriteTObject(outTree, "posteriors", "kOverwrite"); + inFile->Close(); + delete inFile; + + MACH3LOG_INFO("Averaged TTree saved and overwrote original in: {}", TempFilePath); +} From 04f7bb436e75d65ff5385032ac9b58e3104784ac Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Sun, 10 Nov 2024 00:44:03 +0000 Subject: [PATCH 2/4] some tidy --- Diagnostics/ProcessMCMC.cpp | 91 +++++++++++++------------------ mcmc/MCMCProcessor.cpp | 104 ++++++++++++------------------------ mcmc/MCMCProcessor.h | 9 ++-- 3 files changed, 77 insertions(+), 127 deletions(-) diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 37b3b34df..57c839a7f 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -2,21 +2,26 @@ #include "mcmc/MCMCProcessor.h" #include "manager/manager.h" +/// @file ProcessMCMC.cpp +/// @brief Main exectable responsible for different types of MCMC processing like drawing posteriors, triangle plots etc. Actual implantation of methods is in MCMCProcessor + /// @brief Main function processing MCMC and Producing plots inline void ProcessMCMC(const std::string& inputFile); /// @brief Function producing comparison of posterior and more betwen a few MCMC chains inline void MultipleProcessMCMC(); -inline void CalcBayesFactor(MCMCProcessor* Processor); -inline void CalcSavageDickey(MCMCProcessor* Processor); -inline void CalcBipolarPlot(MCMCProcessor* Processor); -inline void CalcParameterEvolution(MCMCProcessor* Processor); -inline void GetTrianglePlot(MCMCProcessor* Processor); -inline void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile); -inline void ReweightPrior(MCMCProcessor* Processor); +inline void CalcBayesFactor(const std::unique_ptr& Processor); +inline void CalcSavageDickey(const std::unique_ptr& Processor); +inline void CalcBipolarPlot(const std::unique_ptr& Processor); +inline void CalcParameterEvolution(const std::unique_ptr& Processor); +inline void GetTrianglePlot(const std::unique_ptr& Processor); +inline void DiagnoseCovarianceMatrix(const std::unique_ptr& Processor, const std::string& inputFile); +inline void ReweightPrior(const std::unique_ptr& Processor); /// @brief KS: Convert TMatrix to TH2D, mostly useful for making fancy plots inline TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title); /// @brief KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution -inline void KolmogorovSmirnovTest(const std::vector>& Processor, TCanvas* Posterior, TString canvasname); +inline void KolmogorovSmirnovTest(const std::vector>& Processor, + const std::unique_ptr& Posterior, + const TString& canvasname); int nFiles; std::vector FileNames; @@ -60,7 +65,6 @@ int main(int argc, char *argv[]) MultipleProcessMCMC(); } - return 0; } @@ -68,12 +72,12 @@ void ProcessMCMC(const std::string& inputFile) { MACH3LOG_INFO("File for study: {} with config {}", inputFile, config); // Make the processor) - MCMCProcessor* Processor = new MCMCProcessor(inputFile); + auto Processor = std::make_unique(inputFile); YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; - bool PlotCorr = GetFromManager(Settings["PlotCorr"], false); + const bool PlotCorr = GetFromManager(Settings["PlotCorr"], false); Processor->SetExcludedTypes(GetFromManager>(Settings["ExcludedTypes"], {""})); Processor->SetExcludedNames(GetFromManager>(Settings["ExcludedNames"], {""})); @@ -145,8 +149,6 @@ void ProcessMCMC(const std::string& inputFile) if(GetFromManager(Settings["DiagnoseCovarianceMatrix"], false)) DiagnoseCovarianceMatrix(Processor, inputFile); } if(GetFromManager(Settings["ReweightPrior"], false)) ReweightPrior(Processor); - - delete Processor; } void MultipleProcessMCMC() @@ -184,7 +186,7 @@ void MultipleProcessMCMC() } // Open a TCanvas to write the posterior onto - TCanvas* Posterior = new TCanvas("PosteriorMulti", "PosteriorMulti", 0, 0, 1024, 1024); + auto Posterior = std::make_unique("PosteriorMulti", "PosteriorMulti", 0, 0, 1024, 1024); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); Posterior->SetGrid(); @@ -214,7 +216,7 @@ void MultipleProcessMCMC() { // This holds the posterior density std::vector hpost(nFiles); - std::vector hpd(nFiles); + std::vector> hpd(nFiles); hpost[0] = static_cast(Processor[0]->GetHpost(i)->Clone()); bool Skip = false; @@ -277,7 +279,8 @@ void MultipleProcessMCMC() TString rebinLeg = Form("#splitline{%s}{#mu = %.2f, #sigma = %.2f}", TitleNames[ik].c_str(), hpost[ik]->GetMean(), hpost[ik]->GetRMS()); leg->AddEntry(hpost[ik], rebinLeg, "l"); - hpd[ik] = new TLine(hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMinimum(), hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMaximum()); + hpd[ik] = std::make_unique(hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMinimum(), + hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMaximum()); hpd[ik]->SetLineColor(hpost[ik]->GetLineColor()); hpd[ik]->SetLineWidth(2); hpd[ik]->SetLineStyle(kSolid); @@ -295,10 +298,8 @@ void MultipleProcessMCMC() leg->Draw("same"); Posterior->cd(); Posterior->Print(canvasname); - for (int ik = 0; ik < nFiles; ik++) - { + for (int ik = 0; ik < nFiles; ik++) { delete hpost[ik]; - delete hpd[ik]; } }//End loop over parameters @@ -313,12 +314,10 @@ void MultipleProcessMCMC() MACH3LOG_INFO("Closing pdf {}", canvasname); canvasname+="]"; Posterior->Print(canvasname); - - delete Posterior; } // KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials -void CalcBayesFactor(MCMCProcessor* Processor) +void CalcBayesFactor(const std::unique_ptr& Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; @@ -336,10 +335,9 @@ void CalcBayesFactor(MCMCProcessor* Processor) } Processor->GetBayesFactor(ParNames, Model1Bounds, Model2Bounds, ModelNames); - return; } -void CalcSavageDickey(MCMCProcessor* Processor) +void CalcSavageDickey(const std::unique_ptr& Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; @@ -354,45 +352,38 @@ void CalcSavageDickey(MCMCProcessor* Processor) EvaluationPoint.push_back(d[1].as()); Bounds.push_back(d[2].as>()); } - Processor->GetSavageDickey(ParNames, EvaluationPoint, Bounds); - return; } -void CalcParameterEvolution(MCMCProcessor* Processor) +void CalcParameterEvolution(const std::unique_ptr& Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; std::vector Intervals; - for (const auto& d : Settings["ParameterEvolution"]) { ParNames.push_back(d[0].as()); Intervals.push_back(d[1].as()); } Processor->ParameterEvolution(ParNames, Intervals); - return; } -void CalcBipolarPlot(MCMCProcessor* Processor) +void CalcBipolarPlot(const std::unique_ptr& Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; - for (const auto& d : Settings["BipolarPlot"]) { ParNames.push_back(d[0].as()); } Processor->GetPolarPlot(ParNames); - return; } - -void GetTrianglePlot(MCMCProcessor* Processor) { +void GetTrianglePlot(const std::unique_ptr& Processor) { YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; @@ -412,12 +403,12 @@ void GetTrianglePlot(MCMCProcessor* Processor) { } /// @brief KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable -void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile) +void DiagnoseCovarianceMatrix(const std::unique_ptr& Processor, const std::string& inputFile) { //Turn of plots from Processor Processor->SetPrintToPDF(false); // Open a TCanvas to write the posterior onto - TCanvas* Canvas = new TCanvas("Canvas", "Canvas", 0, 0, 1024, 1024); + auto Canvas = std::make_unique("Canvas", "Canvas", 0, 0, 1024, 1024); Canvas->SetGrid(); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); @@ -585,17 +576,14 @@ void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& input if(CorrelationPreviousHist != nullptr) delete CorrelationPreviousHist; if(CovarianceHist != nullptr) delete CovarianceHist; if(CorrelationHist != nullptr) delete CorrelationHist; - delete Canvas; } -void ReweightPrior(MCMCProcessor* Processor) +void ReweightPrior(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; const auto& Prior = Settings["PriorReweighting"]; - std::vector Names = Prior[0].as>(); std::vector NewCentral = Prior[1].as>(); std::vector NewError = Prior[2].as>(); @@ -619,7 +607,9 @@ TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title) } //KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution -void KolmogorovSmirnovTest(const std::vector>& Processor, TCanvas* Posterior, TString canvasname) +void KolmogorovSmirnovTest(const std::vector>& Processor, + const std::unique_ptr& Posterior, + const TString& canvasname) { const Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2}; const Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted}; @@ -627,8 +617,8 @@ void KolmogorovSmirnovTest(const std::vector>& Pr for(int i = 0; i < Processor[0]->GetNParams(); ++i) { // This holds the posterior density - TH1D **hpost = new TH1D*[nFiles]; - TH1D **CumulativeDistribution = new TH1D*[nFiles]; + std::vector hpost(nFiles); + std::vector CumulativeDistribution(nFiles); TString Title; double Prior = 1.0; @@ -675,8 +665,6 @@ void KolmogorovSmirnovTest(const std::vector>& Pr delete hpost[ik]; delete CumulativeDistribution[ik]; } - delete[] hpost; - delete[] CumulativeDistribution; continue; } @@ -696,15 +684,15 @@ void KolmogorovSmirnovTest(const std::vector>& Pr std::vector TestStatBin(nFiles, 0); std::vector TestStatD(nFiles, -999); - std::vector LineD(nFiles); + std::vector> LineD(nFiles); //Find KS statistic for (int ik = 1 ; ik < nFiles; ik++) { const int NumberOfBins = CumulativeDistribution[0]->GetXaxis()->GetNbins(); for (int j = 1; j < NumberOfBins+1; ++j) { - double BinValue = CumulativeDistribution[0]->GetBinCenter(j); - int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue); + const double BinValue = CumulativeDistribution[0]->GetBinCenter(j); + const int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue); //KS: Calculate D statistic for this bin, only save it if it's bigger than previously found value double TempDstat = std::fabs(CumulativeDistribution[0]->GetBinContent(j) - CumulativeDistribution[ik]->GetBinContent(BinNumber)); if(TempDstat > TestStatD[ik]) @@ -717,7 +705,7 @@ void KolmogorovSmirnovTest(const std::vector>& Pr for (int ik = 0 ; ik < nFiles; ik++) { - LineD[ik] = new TLine(CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), 0, CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), CumulativeDistribution[0]->GetBinContent(TestStatBin[ik])); + LineD[ik] = std::make_unique(CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), 0, CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), CumulativeDistribution[0]->GetBinContent(TestStatBin[ik])); LineD[ik]->SetLineColor(CumulativeColor[ik]); LineD[ik]->SetLineWidth(2.0); } @@ -730,7 +718,7 @@ void KolmogorovSmirnovTest(const std::vector>& Pr for (int ik = 0; ik < nFiles; ik++) leg->AddEntry(CumulativeDistribution[ik], TitleNames[ik].c_str(), "l"); for (int ik = 1; ik < nFiles; ik++) - leg->AddEntry(LineD[ik], Form("#Delta D = %.4f", TestStatD[ik]), "l"); + leg->AddEntry(LineD[ik].get(), Form("#Delta D = %.4f", TestStatD[ik]), "l"); leg->SetLineColor(0); leg->SetLineStyle(0); @@ -748,9 +736,6 @@ void KolmogorovSmirnovTest(const std::vector>& Pr { delete hpost[ik]; delete CumulativeDistribution[ik]; - delete LineD[ik]; } - delete[] hpost; - delete[] CumulativeDistribution; } //End loop over parameter } diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index 46f14a73a..130c1cfbe 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -25,7 +25,6 @@ MCMCProcessor::MCMCProcessor(const std::string &InputFile) : StepNumber = nullptr; Posterior = nullptr; - hpost = nullptr; hpost2D = nullptr; hviolin = nullptr; hviolin_prior = nullptr; @@ -121,13 +120,10 @@ MCMCProcessor::~MCMCProcessor() { delete Errors_HPD_Positive; delete Errors_HPD_Negative; - if(hpost != nullptr) + + for (int i = 0; i < nDraw; ++i) { - for (int i = 0; i < nDraw; ++i) - { - delete hpost[i]; - } - delete[] hpost; + if(hpost[i] != nullptr) delete hpost[i]; } if(CacheMCMC) { @@ -1202,7 +1198,6 @@ void MCMCProcessor::CacheSteps() { MACH3LOG_INFO("Caching steps took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries ); } - // ********************* // Make the post-fit covariance matrix in all dimensions void MCMCProcessor::MakeCovariance_MP(bool Mute) { @@ -1309,7 +1304,6 @@ void MCMCProcessor::MakeCovariance_MP(bool Mute) { } } - // ********************* // Based on @cite roberts2009adaptive // all credits for finding and studying it goes to Henry @@ -2298,8 +2292,7 @@ void MCMCProcessor::SetupOutput() { (*Correlation)(i, j) = _UNDEF_; } } - - hpost = new TH1D*[nDraw](); + hpost.resize(nDraw); } // **************************** @@ -2371,7 +2364,7 @@ TH1D* MCMCProcessor::MakePrefit() { PreFitPlot->SetBinError(i+1, Error); PreFitPlot->GetXaxis()->SetBinLabel(i+1, ParamNames[ParamEnum][ParamNo]); } - PreFitPlot->SetDirectory(0); + PreFitPlot->SetDirectory(nullptr); PreFitPlot->SetFillStyle(1001); PreFitPlot->SetFillColor(kRed-3); @@ -2623,7 +2616,6 @@ void MCMCProcessor::ReadFDFile() { // Read the Osc cov file and get the input central values and errors void MCMCProcessor::ReadOSCFile() { // *************** - YAML::Node OscFile = CovConfig[kOSCPar];; auto systematics = OscFile["Systematics"]; @@ -2893,7 +2885,6 @@ void MCMCProcessor::GetSavageDickey(const std::vector& ParNames, const std::vector& EvaluationPoint, const std::vector>& Bounds){ // ************************** - if((ParNames.size() != EvaluationPoint.size()) || (Bounds.size() != EvaluationPoint.size())) { MACH3LOG_ERROR("Size doesn't match"); @@ -3037,7 +3028,6 @@ void MCMCProcessor::ReweightPrior(const std::vector& Names, const std::vector& NewCentral, const std::vector& NewError) { // ************************** - MACH3LOG_INFO("Reweighting Prior"); if( (Names.size() != NewCentral.size()) || (NewCentral.size() != NewError.size())) @@ -3167,7 +3157,8 @@ void MCMCProcessor::ParameterEvolution(const std::vector& Names, for(int i = NIntervals[k]-1; i >= 0; --i) { // This holds the posterior density - TH1D* EvePlot = new TH1D(BranchNames[ParamNo], BranchNames[ParamNo], nBins, mini, maxi); + auto EvePlot = std::make_unique(BranchNames[ParamNo], BranchNames[ParamNo], nBins, mini, maxi); + EvePlot->SetDirectory(nullptr); EvePlot->SetMinimum(0); EvePlot->GetYaxis()->SetTitle("PDF"); EvePlot->GetYaxis()->SetNoExponent(false); @@ -3190,17 +3181,14 @@ void MCMCProcessor::ParameterEvolution(const std::vector& Names, EvePlot->Draw("HIST"); - TText *text = new TText(0.3, 0.8, TextTitle.c_str()); - text->SetTextFont (43); - text->SetTextSize (40); - text->SetNDC(true); - text->Draw("SAME"); + TText text(0.3, 0.8, TextTitle.c_str()); + text.SetTextFont (43); + text.SetTextSize (40); + text.SetNDC(true); + text.Draw("SAME"); if(i == 0) Posterior->Print(((Names[k] + ".gif++20").c_str())); // produces infinite loop animated GIF else Posterior->Print(((Names[k]+".gif+20").c_str())); // add picture to .gif - - delete EvePlot; - delete text; Counter++; } } @@ -3405,13 +3393,12 @@ void MCMCProcessor::PrepareDiagMCMC() { //CW: Draw trace plots of the parameters i.e. parameter vs step void MCMCProcessor::ParamTraces() { // ***************** - if (ParStep == nullptr) PrepareDiagMCMC(); MACH3LOG_INFO("Making trace plots..."); // Make the TH1Ds - TH1D** TraceParamPlots = new TH1D*[nDraw]; - TH1D** TraceSamplePlots = new TH1D*[nSamples]; - TH1D** TraceSystsPlots = new TH1D*[nSysts]; + std::vector TraceParamPlots(nDraw); + std::vector TraceSamplePlots(nSamples); + std::vector TraceSystsPlots(nSysts); // Set the titles and limits for TH2Ds for (int j = 0; j < nDraw; ++j) { @@ -3422,7 +3409,6 @@ void MCMCProcessor::ParamTraces() { GetNthParameter(j, Prior, PriorError, Title); std::string HistName = Form("%s_%s_Trace", Title.Data(), BranchNames[j].Data()); - TraceParamPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nEntries, 0, nEntries); TraceParamPlots[j]->GetXaxis()->SetTitle("Step"); TraceParamPlots[j]->GetYaxis()->SetTitle("Parameter Variation"); @@ -3443,7 +3429,6 @@ void MCMCProcessor::ParamTraces() { } // Have now made the empty TH1Ds, now for writing content to them! - // Loop over the number of parameters to draw their traces // Each histogram #ifdef MULTITHREAD @@ -3455,11 +3440,9 @@ void MCMCProcessor::ParamTraces() { for (int j = 0; j < nDraw; ++j) { TraceParamPlots[j]->SetBinContent(i, ParStep[j][i]); } - for (int j = 0; j < nSamples; ++j) { TraceSamplePlots[j]->SetBinContent(i, SampleValues[i][j]); } - for (int j = 0; j < nSysts; ++j) { TraceSystsPlots[j]->SetBinContent(i, SystValues[i][j]); } @@ -3477,7 +3460,6 @@ void MCMCProcessor::ParamTraces() { delete Fitter; delete TraceParamPlots[j]; } - delete[] TraceParamPlots; TDirectory *LLDir = OutputFile->mkdir("LogL"); LLDir->cd(); @@ -3486,7 +3468,6 @@ void MCMCProcessor::ParamTraces() { delete TraceSamplePlots[j]; delete[] SampleValues[j]; } - delete[] TraceSamplePlots; delete[] SampleValues; for (int j = 0; j < nSysts; ++j) { @@ -3494,7 +3475,6 @@ void MCMCProcessor::ParamTraces() { delete TraceSystsPlots[j]; delete SystValues[j]; } - delete[] TraceSystsPlots; delete[] SystValues; TraceDir->Close(); @@ -3507,7 +3487,6 @@ void MCMCProcessor::ParamTraces() { //KS: Calculate autocorrelations supports both OpenMP and CUDA :) void MCMCProcessor::AutoCorrelation() { // ********************************* - if (ParStep == nullptr) PrepareDiagMCMC(); TStopwatch clock; @@ -3516,15 +3495,15 @@ void MCMCProcessor::AutoCorrelation() { MACH3LOG_INFO("Making auto-correlations for nLags = {}", nLags); // The sum of (Y-Ymean)^2 over all steps for each parameter - double **DenomSum = new double*[nDraw](); - double **NumeratorSum = new double*[nDraw](); - double **LagL = new double*[nDraw]; + std::vector> DenomSum(nDraw); + std::vector> NumeratorSum(nDraw); + std::vector> LagL(nDraw); for (int j = 0; j < nDraw; ++j) { - DenomSum[j] = new double[nLags]; - NumeratorSum[j] = new double[nLags]; - LagL[j] = new double[nLags]; + DenomSum[j].resize(nLags); + NumeratorSum[j].resize(nLags); + LagL[j].resize(nLags); } - TH1D** LagKPlots = new TH1D*[nDraw]; + std::vector LagKPlots(nDraw); // Loop over the parameters of interest for (int j = 0; j < nDraw; ++j) { @@ -3566,7 +3545,6 @@ void MCMCProcessor::AutoCorrelation() { const double Product = Diff*LagTerm; NumeratorSum[j][k] += Product; } - // Square the difference to form the denominator const double Denom = Diff*Diff; DenomSum[j][k] += Denom; @@ -3628,19 +3606,10 @@ void MCMCProcessor::AutoCorrelation() { LagKPlots[j]->Write(); delete LagKPlots[j]; } - delete[] LagKPlots; //KS: This is different diagnostic however it relies on calculated Lag, thus we call it before we delete LagKPlots CalculateESS(nLags, LagL); - for (int j = 0; j < nDraw; ++j) { - delete[] NumeratorSum[j]; - delete[] DenomSum[j]; - delete[] LagL[j]; - } - delete[] NumeratorSum; - delete[] DenomSum; - delete[] LagL; delete[] ParamSums; AutoCorrDir->Close(); @@ -3737,11 +3706,11 @@ void MCMCProcessor::PrepareGPU_AutoCorr(const int nLags) { // KS: calc Effective Sample Size Following @cite StanManual // Furthermore we calculate Sampling efficiency following @cite hanson2008mcmc // Rule of thumb is to have efficiency above 25% -void MCMCProcessor::CalculateESS(const int nLags, double** LagL) { +void MCMCProcessor::CalculateESS(const int nLags, const std::vector>& LagL) { // ************************** - if(LagL == nullptr) + if(LagL.size() == 0) { - MACH3LOG_ERROR("LagL is nullptr"); + MACH3LOG_ERROR("Size of LagL is 0"); throw MaCh3Exception(__FILE__ , __LINE__ ); } MACH3LOG_INFO("Making ESS plots..."); @@ -3842,11 +3811,10 @@ void MCMCProcessor::CalculateESS(const int nLags, double** LagL) { //CW: Batched means, literally read from an array and chuck into TH1D void MCMCProcessor::BatchedMeans() { // ************************** - if (BatchedAverages == nullptr) PrepareDiagMCMC(); MACH3LOG_INFO("Making BatchedMeans plots..."); - TH1D ** BatchedParamPlots = new TH1D*[nDraw]; + std::vector BatchedParamPlots(nDraw); for (int j = 0; j < nDraw; ++j) { TString Title = ""; double Prior = 1.0; @@ -3882,7 +3850,6 @@ void MCMCProcessor::BatchedMeans() { delete Fitter; delete BatchedParamPlots[j]; } - delete[] BatchedParamPlots; //KS: Get the batched means variance estimation and variable indicating if number of batches is sensible // We do this before deleting BatchedAverages @@ -4077,11 +4044,10 @@ void MCMCProcessor::PowerSpectrumAnalysis() { TDirectory *PowerDir = OutputFile->mkdir("PowerSpectrum"); PowerDir->cd(); - TGraph **plot = new TGraph*[nPrams]; TVectorD* PowerSpectrumStepSize = new TVectorD(nPrams); for (int j = 0; j < nPrams; ++j) { - plot[j] = new TGraph(v_size, k_j[j].data(), P_j[j].data()); + TGraph* plot = new TGraph(v_size, k_j[j].data(), P_j[j].data()); TString Title = ""; double Prior = 1.0; @@ -4091,10 +4057,10 @@ void MCMCProcessor::PowerSpectrumAnalysis() { std::string name = Form("Power Spectrum of %s;k;P(k)", Title.Data()); - plot[j]->SetTitle(name.c_str()); + plot->SetTitle(name.c_str()); name = Form("%s_power_spectrum", Title.Data()); - plot[j]->SetName(name.c_str()); - plot[j]->SetMarkerStyle(7); + plot->SetName(name.c_str()); + plot->SetMarkerStyle(7); // Equation 18 TF1 *func = new TF1("power_template", "[0]*( ([1] / x)^[2] / (([1] / x)^[2] +1) )", 0.0, 1.0); @@ -4110,23 +4076,21 @@ void MCMCProcessor::PowerSpectrumAnalysis() { func->SetParLimits(1, 0.001, 1.0); // k* should be within a reasonable range func->SetParLimits(2, 0.0, 5.0); // alpha should be positive - plot[j]->Fit("power_template","Rq"); + plot->Fit("power_template","Rq"); Posterior->SetLogx(); Posterior->SetLogy(); Posterior->SetGrid(); - plot[j]->Write(plot[j]->GetName()); - plot[j]->Draw("AL"); + plot->Write(plot->GetName()); + plot->Draw("AL"); func->Draw("SAME"); if(printToPDF) Posterior->Print(CanvasName); //KS: I have no clue what is the reason behind this. Found this in Rick Calland code... (*PowerSpectrumStepSize)(j) = std::sqrt(func->GetParameter(0)/float(v_size*0.5)); - delete func; - delete plot[j]; + delete plot; } - delete [] plot; PowerSpectrumStepSize->Write("PowerSpectrumStepSize"); delete PowerSpectrumStepSize; diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 6b6143357..7be5086c0 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -253,6 +253,7 @@ class MCMCProcessor { /// @brief Sett output suffix, this way jobs using the same file will have different names inline void SetOutputSuffix(const std::string Suffix){OutputSuffix = Suffix; }; + /// @brief Allow to set addtional cuts based on ROOT TBrowser cut, for to only affect one mass ordering inline void SetPosterior1DCut(const std::string Cut){Posterior1DCut = Cut; }; private: /// @brief Prepare prefit histogram for parameter overlay plot @@ -302,7 +303,7 @@ class MCMCProcessor { /// \cite StanManual /// \cite hanson2008mcmc /// \cite gabry2024visual - inline void CalculateESS(const int nLags, double **LagL); + inline void CalculateESS(const int nLags, const std::vector>& LagL); /// @brief Get the batched means variance estimation and variable indicating if number of batches is sensible /// \cite chakraborty2019estimating /// \cite rossetti2024batch @@ -315,7 +316,8 @@ class MCMCProcessor { inline void GewekeDiagnostic(); /// @brief Acceptance Probability inline void AcceptanceProbabilities(); - /// @brief RC: Perform spectral analysis of MCMC based on \cite Dunkley:2004sv + /// @brief RC: Perform spectral analysis of MCMC + /// \cite Dunkley:2004sv inline void PowerSpectrumAnalysis(); /// Name of MCMC file @@ -327,7 +329,6 @@ class MCMCProcessor { /// Covariance matrix config std::vector CovConfig; - /// Main chain storing all steps etc TChain *Chain; /// BurnIn Cuts @@ -449,7 +450,7 @@ class MCMCProcessor { TMatrixDSym *Correlation; /// Holds 1D Posterior Distributions - TH1D **hpost; + std::vector hpost; /// Holds 2D Posterior Distributions TH2D ***hpost2D; /// Holds violin plot for all dials From 258501679649daf6098f221d389597b428f3c414 Mon Sep 17 00:00:00 2001 From: Kamil <45295406+KSkwarczynski@users.noreply.github.com> Date: Thu, 14 Nov 2024 00:21:22 +0000 Subject: [PATCH 3/4] Add comment about potential impact on tail of distuvtion --- mcmc/StatisticalUtils.h | 1 + 1 file changed, 1 insertion(+) diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h index 94a9ae225..510eff886 100644 --- a/mcmc/StatisticalUtils.h +++ b/mcmc/StatisticalUtils.h @@ -593,6 +593,7 @@ inline void ThinningMCMC(const std::string& FilePath, const int ThinningCut) { /// @param AveragingCut Number of consecutive entries to average /// @cite nesterov2009primal /// @warning Averaging is done over entries, so this may not work well for merged chains or non-sequential entries. +/// @warning Impact on tails in very long chains haven't been checked, thus be carefull /// @note The resulting file overwrites the original file. inline void AveragingMCMC(const std::string& FilePath, const int AveragingCut) { // ******************** From b79c9d344699c39b1aea17ecd20f9461fc94cdc6 Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Sat, 16 Nov 2024 00:55:22 +0000 Subject: [PATCH 4/4] remove averaging --- Diagnostics/ProcessMCMC.cpp | 5 +-- Doc/bibliography.bib | 11 ----- mcmc/MCMCProcessor.h | 8 +--- mcmc/StatisticalUtils.h | 90 ------------------------------------- 4 files changed, 3 insertions(+), 111 deletions(-) diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 57c839a7f..7abca6cbf 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -99,9 +99,8 @@ void ProcessMCMC(const std::string& inputFile) Processor->Initialise(); if(Settings["Thinning"]) { - bool DoThin = Settings["Thinning"][0].as(); - if(DoThin){ - Processor->ThinMCMC(Settings["Thinning"][1].as(), GetFromManager(Settings["Average"], false)); + if(Settings["Thinning"][0].as()){ + Processor->ThinMCMC(Settings["Thinning"][1].as()); } } // Make the postfit diff --git a/Doc/bibliography.bib b/Doc/bibliography.bib index 35ec3846e..79fac0844 100755 --- a/Doc/bibliography.bib +++ b/Doc/bibliography.bib @@ -180,14 +180,3 @@ @article{2011ThinningMCMC doi = {10.1111/j.2041-210X.2011.00131.x}, url = {https://doi.org/10.1111/j.2041-210X.2011.00131.x}, } - -@article{nesterov2009primal, - author = {Nesterov, Yurii}, - title = {Primal-dual subgradient methods for convex problems}, - journal = {Mathematical Programming}, - volume = {120}, - pages = {221--259}, - year = {2009}, - doi = {10.1007/s10107-007-0149-x}, - note = {Received 29 September 2005; Accepted 13 January 2007; Published 19 June 2007; Issue Date August 2009} -} diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 7be5086c0..ef91d7a00 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -162,13 +162,7 @@ class MCMCProcessor { /// @brief Thin MCMC Chain, to save space and maintain low autocorrelations. /// @param ThinningCut every which entry you want to thin /// @param Average If true will perform MCMC averaging instead of thinning - inline void ThinMCMC(const int ThinningCut, const bool Average = false) { - if(Average){ - AveragingMCMC(MCMCFile+".root", ThinningCut); - } else { - ThinningMCMC(MCMCFile+".root", ThinningCut); - } - }; + inline void ThinMCMC(const int ThinningCut) { ThinningMCMC(MCMCFile+".root", ThinningCut); }; /// @brief KS: Perform MCMC diagnostic including Autocorrelation, Trace etc. void DiagMCMC(); diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h index 510eff886..4312b6b3a 100644 --- a/mcmc/StatisticalUtils.h +++ b/mcmc/StatisticalUtils.h @@ -584,93 +584,3 @@ inline void ThinningMCMC(const std::string& FilePath, const int ThinningCut) { MACH3LOG_INFO("Thinned TTree saved and overwrote original in: {}", TempFilePath); } - - -// ******************** -/// @brief Average MCMC Chain entries, to reduce autocorrelation and improve estimate stability. -/// -/// @param FilePath Path to the MCMC chain you want to average -/// @param AveragingCut Number of consecutive entries to average -/// @cite nesterov2009primal -/// @warning Averaging is done over entries, so this may not work well for merged chains or non-sequential entries. -/// @warning Impact on tails in very long chains haven't been checked, thus be carefull -/// @note The resulting file overwrites the original file. -inline void AveragingMCMC(const std::string& FilePath, const int AveragingCut) { -// ******************** - if (AveragingCut <= 1) { - MACH3LOG_WARN("Invalid AveragingCut value ({}). It must be greater than 1.", AveragingCut); - throw MaCh3Exception(__FILE__, __LINE__); - } - - // Define the path for the temporary averaged file - std::string TempFilePath = "Averaged_" + FilePath; - int ret = system(("cp " + FilePath + " " + TempFilePath).c_str()); - if (ret != 0) { - MACH3LOG_WARN("Error: system call to copy file failed with code {}", ret); - } - - TFile *inFile = TFile::Open(TempFilePath.c_str(), "UPDATE"); - if (!inFile || inFile->IsZombie()) { - MACH3LOG_ERROR("Error opening file: {}", TempFilePath); - throw MaCh3Exception(__FILE__, __LINE__); - } - - TTree *inTree = inFile->Get("posteriors"); - if (!inTree) { - MACH3LOG_ERROR("Error: TTree 'posteriors' not found in file."); - inFile->ls(); - inFile->Close(); - throw MaCh3Exception(__FILE__, __LINE__); - } - - // Clone the structure without data - TTree *outTree = inTree->CloneTree(0); - Long64_t nEntries = inTree->GetEntries(); - int nBranches = inTree->GetNbranches(); - - std::vector AvgValues(nBranches, 0.0); - std::vector BranchData(nBranches, 0.0); - - // Loop over branches and set up data pointers - TObjArray *branches = inTree->GetListOfBranches(); - for (int b = 0; b < branches->GetEntries(); ++b) { - TBranch *branch = static_cast(branches->At(b)); - inTree->SetBranchAddress(branch->GetName(), &BranchData[b]); - } - - MACH3LOG_INFO("Averaging every {} entries. Total entries: {}", AveragingCut, nEntries); - - int batchCount = 0; - for (Long64_t i = 0; i < nEntries; i += AveragingCut) { - // Reset to 0 - std::fill(AvgValues.begin(), AvgValues.end(), 0.0); - - /// @todo Here we perform normal mean. But in literature there are available weighted means and more complex operations - int entriesInBatch = 0; - for (int j = 0; j < AveragingCut && (i + j) < nEntries; ++j) { - inTree->GetEntry(i + j); - for (int b = 0; b < branches->GetEntries(); ++b) { - AvgValues[b] += BranchData[b]; - } - ++entriesInBatch; - } - for (int b = 0; b < branches->GetEntries(); ++b) { - AvgValues[b] /= entriesInBatch; - BranchData[b] = AvgValues[b]; - } - - outTree->Fill(); - - if (i % (nEntries/10) == 0) { - MaCh3Utils::PrintProgressBar(batchCount, nEntries / AveragingCut); - } - ++batchCount; - } - - // Write the averaged tree back to the file - inFile->WriteTObject(outTree, "posteriors", "kOverwrite"); - inFile->Close(); - delete inFile; - - MACH3LOG_INFO("Averaged TTree saved and overwrote original in: {}", TempFilePath); -}