From 395c63cc0a125646fbbfbe3917f5b166dbf011f7 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Tue, 19 Nov 2024 17:27:24 -0500 Subject: [PATCH 01/13] Implemented a faster / lower memory RHat calculator and moved the standard RHat into a "HighMem" macro --- Diagnostics/RHat.cpp | 332 ++++---------- Diagnostics/RHat_HighMem.cpp | 845 +++++++++++++++++++++++++++++++++++ 2 files changed, 934 insertions(+), 243 deletions(-) create mode 100644 Diagnostics/RHat_HighMem.cpp diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 86ffc76d0..3be941df5 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -29,18 +29,31 @@ /// and helps determine whether the chains have reached a stable distribution. /// /// @cite gelman2019. +/// +/// MJR: Update -- Improved memory usage so that whole chains can be quickly loaded without requiring copious amounts +/// of RAM. This comes at the cost of not being able to calculate the Folded RHat since finding the median +/// requires the loading of full chains at a time. The method has been validated to give identical results to the +/// "High Memory" (original) version at a fraction of the runtime and resources. +/// +/// The input format is also slightly altered; since we can now load entire chains, there's less need to +/// specify how many toys are desired for a sub-sample, so the Ntoys input has been removed. +/// // ******************* -int Ntoys; +int* Ntoys; +int TotToys; +int NThin; int Nchains; - int nDraw; std::vector BranchNames; std::vector MCMCFile; std::vector ValidPar; -double ***Draws; +double* S1_global; // Sum_i^N x_i | total +double* S2_global; // Sum_i^N x_i^2 | total +double** S1_chain; // Sum_i^N x_i | for each chain +double** S2_chain; // Sum_i^N x_i^2 | for each chain double** Mean; double** StandardDeviation; @@ -53,19 +66,6 @@ double* MarginalPosteriorVariance; double* RHat; double* EffectiveSampleSize; -double ***DrawsFolded; -double* MedianArr; - -double** MeanFolded; -double** StandardDeviationFolded; - -double* MeanGlobalFolded; -double* StandardDeviationGlobalFolded; - -double* BetweenChainVarianceFolded; -double* MarginalPosteriorVarianceFolded; -double* RHatFolded; -double* EffectiveSampleSizeFolded; // ******************* void PrepareChains(); void InitialiseArrays(); @@ -86,7 +86,6 @@ int main(int argc, char *argv[]) { SetMaCh3LoggerFormat(); MaCh3Utils::MaCh3Welcome(); - Draws = nullptr; Mean = nullptr; StandardDeviation = nullptr; @@ -98,50 +97,31 @@ int main(int argc, char *argv[]) { RHat = nullptr; EffectiveSampleSize = nullptr; - DrawsFolded = nullptr; - MedianArr = nullptr; - MeanFolded = nullptr; - StandardDeviationFolded = nullptr; - - MeanGlobalFolded = nullptr; - StandardDeviationGlobalFolded = nullptr; - - BetweenChainVarianceFolded = nullptr; - MarginalPosteriorVarianceFolded = nullptr; - RHatFolded = nullptr; - EffectiveSampleSizeFolded = nullptr; - Nchains = 0; - if (argc == 1 || argc == 2) + if (argc < 2) { MACH3LOG_ERROR("Wrong arguments"); - MACH3LOG_ERROR("./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); + MACH3LOG_ERROR("./RHat MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); throw MaCh3Exception(__FILE__ , __LINE__ ); } Ntoys = atoi(argv[1]); //KS Gelman suggests to diagnose on more than one chain - for (int i = 2; i < argc; i++) + for (int i = 1; i < argc; i++) { MCMCFile.push_back(std::string(argv[i])); MACH3LOG_INFO("Adding file: {}", MCMCFile.back()); Nchains++; } - if(Ntoys < 1) - { - MACH3LOG_ERROR("You specified {} specify larger greater than 0", Ntoys); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - if(Nchains == 1) { MACH3LOG_WARN("Gelman is going to be sad :(. He suggested you should use more than one chain (at least 4). Code works fine for one chain, however, estimator might be biased."); MACH3LOG_WARN("Multiple chains are more likely to reveal multimodality and poor adaptation or mixing:"); } - MACH3LOG_INFO("Diagnosing {} chains, with {} toys", Nchains, Ntoys); + MACH3LOG_INFO("Diagnosing {} chains", Nchains); PrepareChains(); @@ -161,20 +141,19 @@ int main(int argc, char *argv[]) { // Load chain and prepare toys void PrepareChains() { // ******************* - auto rnd = std::make_unique(0); - - MACH3LOG_INFO("Generating {}", Ntoys); TStopwatch clock; clock.Start(); + Ntoys = new int[Nchains](); + TotToys = 0; std::vector BurnIn(Nchains); std::vector nEntries(Nchains); std::vector nBranches(Nchains); std::vector step(Nchains); - Draws = new double**[Nchains](); - DrawsFolded = new double**[Nchains](); + S1_chain = new double*[Nchains](); + S2_chain = new double*[Nchains](); // KS: This can reduce time necessary for caching even by half #ifdef MULTITHREAD @@ -187,8 +166,13 @@ void PrepareChains() { { TChain* Chain = new TChain("posteriors"); Chain->Add(MCMCFile[m].c_str()); - MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); + nEntries[m] = int(Chain->GetEntries()); + Ntoys[m] = nEntries[m]; + TotToys += Ntoys[m]; + + MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); + MACH3LOG_INFO("Generating {}", Ntoys); // Set the step cut to be 20% BurnIn[m] = nEntries[m]/5; @@ -249,6 +233,25 @@ void PrepareChains() { if(m == 0) nDraw = int(BranchNames.size()); + // MJR: Initialize quantities needed for calculating RHat + S1_chain[m] = new double[nDraw](); + S2_chain[m] = new double[nDraw](); + if (m == 0) + { + S1_global = new double[nDraw](); + S2_global = new double[nDraw](); + } + for (int id = 0; id < nDraw; ++id) + { + S1_chain[m][id] = 0.0; + S2_chain[m][id] = 0.0; + if (m == 0) + { + S1_global[id] = 0.0; + S2_global[id] = 0.0; + } + } + //TN: Qualitatively faster sanity check, with the very same outcome (all chains have the same #branches) if(m > 0) { @@ -260,18 +263,12 @@ void PrepareChains() { } } - //TN: move the Draws here, so we need to iterate over every chain only once - Draws[m] = new double*[Ntoys](); - DrawsFolded[m] = new double*[Ntoys](); - for(int i = 0; i < Ntoys; i++) + // MJR: Create an array to hold branch values. Resetting branch addresses + // for every step is very expensive. + double* branch_values = new double[nDraw](); + for (int id = 0; id < nDraw; ++id) { - Draws[m][i] = new double[nDraw](); - DrawsFolded[m][i] = new double[nDraw](); - for(int j = 0; j < nDraw; j++) - { - Draws[m][i][j] = 0.; - DrawsFolded[m][i][j] = 0.; - } + Chain->SetBranchAddress(BranchNames[id].Data(), &branch_values[id]); } //TN: move looping over toys here, so we don't need to loop over chains more than once @@ -284,10 +281,11 @@ void PrepareChains() { throw MaCh3Exception(__FILE__ , __LINE__ ); } - for (int i = 0; i < Ntoys; i++) + MACH3LOG_INFO("Loading chain {} / {}...", m, Nchains); + for (int i = 0; i < Ntoys[m]; i++) { - // Get a random entry after burn in - int entry = int(nEntries[m]*rnd->Rndm()); + // This is here as a placeholder in case we want to do some thinning later + int entry = i; Chain->GetEntry(entry); @@ -295,7 +293,7 @@ void PrepareChains() { // Note, entry is not necessarily the same as the step due to merged ROOT files, so can't choose an entry in the range BurnIn - nEntries :( if (step[m] < BurnIn[m]) { - i--; + Ntoys[m]--; continue; } @@ -305,53 +303,25 @@ void PrepareChains() { MACH3LOG_DEBUG("Getting random entry {}", entry); } - // Set the branch addresses for params + // MJR: Fill running quantities instead of loading everything into RAM. + // This is where we save big on both memory and time (resetting + // branch addresses and calling GetEntry() again here is very slow). for (int j = 0; j < nDraw; ++j) { - Chain->SetBranchAddress(BranchNames[j].Data(), &Draws[m][i][j]); + S1_global[j] += branch_values[j]; + S2_global[j] += branch_values[j]*branch_values[j]; + S1_chain[m][j] += branch_values[j]; + S2_chain[m][j] += branch_values[j]*branch_values[j]; } - Chain->GetEntry(entry); }//end loop over toys //TN: There, we now don't need to keep the chain in memory anymore delete Chain; + delete[] branch_values; + MACH3LOG_INFO("Finished loading chain {}!", m); } - //KS: Now prepare folded draws, quoting Gelman - //"We propose to report the maximum of rank normalized split-Rb and rank normalized folded-split-Rb for each parameter" - MedianArr = new double[nDraw](); - #ifdef MULTITHREAD - #pragma omp parallel for - #endif - for(int j = 0; j < nDraw; j++) - { - MedianArr[j] = 0.; - std::vector TempDraws(static_cast(Ntoys) * Nchains); - for(int m = 0; m < Nchains; m++) - { - for(int i = 0; i < Ntoys; i++) - { - const int im = i+m; - TempDraws[im] = Draws[m][i][j]; - } - } - MedianArr[j] = CalcMedian(TempDraws.data(), Ntoys*Nchains); - } - - #ifdef MULTITHREAD - #pragma omp parallel for collapse(3) - #endif - for(int m = 0; m < Nchains; m++) - { - for(int i = 0; i < Ntoys; i++) - { - for(int j = 0; j < nDraw; j++) - { - DrawsFolded[m][i][j] = std::fabs(Draws[m][i][j] - MedianArr[j]); - } - } - } clock.Stop(); MACH3LOG_INFO("Finished calculating Toys, it took {:.2f}s to finish", clock.RealTime()); } @@ -373,31 +343,16 @@ void InitialiseArrays() { RHat = new double[nDraw](); EffectiveSampleSize = new double[nDraw](); - MeanFolded = new double*[Nchains](); - StandardDeviationFolded = new double*[Nchains](); - - MeanGlobalFolded = new double[nDraw](); - StandardDeviationGlobalFolded = new double[nDraw](); - BetweenChainVarianceFolded = new double[nDraw](); - - MarginalPosteriorVarianceFolded = new double[nDraw](); - RHatFolded = new double[nDraw](); - EffectiveSampleSizeFolded = new double[nDraw](); - for (int m = 0; m < Nchains; ++m) { Mean[m] = new double[nDraw](); StandardDeviation[m] = new double[nDraw](); - MeanFolded[m] = new double[nDraw](); - StandardDeviationFolded[m] = new double[nDraw](); for (int j = 0; j < nDraw; ++j) { Mean[m][j] = 0.; StandardDeviation[m][j] = 0.; - MeanFolded[m][j] = 0.; - StandardDeviationFolded[m][j] = 0.; if(m == 0) { MeanGlobal[j] = 0.; @@ -406,13 +361,6 @@ void InitialiseArrays() { MarginalPosteriorVariance[j] = 0.; RHat[j] = 0.; EffectiveSampleSize[j] = 0.; - - MeanGlobalFolded[j] = 0.; - StandardDeviationGlobalFolded[j] = 0.; - BetweenChainVarianceFolded[j] = 0.; - MarginalPosteriorVarianceFolded[j] = 0.; - RHatFolded[j] = 0.; - EffectiveSampleSizeFolded[j] = 0.; } } } @@ -442,7 +390,7 @@ void CalcRhat() { #endif #ifdef MULTITHREAD - #pragma omp for collapse(2) + #pragma omp for #endif //KS: loop over chains and draws are independent so might as well collapse for sweet cache hits //Calculate the mean for each parameter within each considered chain @@ -450,13 +398,8 @@ void CalcRhat() { { for (int j = 0; j < nDraw; ++j) { - for(int i = 0; i < Ntoys; i++) - { - Mean[m][j] += Draws[m][i][j]; - MeanFolded[m][j] += DrawsFolded[m][i][j]; - } - Mean[m][j] = Mean[m][j]/Ntoys; - MeanFolded[m][j] = MeanFolded[m][j]/Ntoys; + Mean[m][j] += S1_chain[m][j] / (double)Ntoys[m]; + StandardDeviation[m][j] = S2_chain[m][j]/(double)Ntoys[m] - Mean[m][j]*Mean[m][j]; } } @@ -468,45 +411,14 @@ void CalcRhat() { { for (int m = 0; m < Nchains; ++m) { - MeanGlobal[j] += Mean[m][j]; - MeanGlobalFolded[j] += MeanFolded[m][j]; - } - MeanGlobal[j] = MeanGlobal[j]/Nchains; - MeanGlobalFolded[j] = MeanGlobalFolded[j]/Nchains; - } - - - #ifdef MULTITHREAD - #pragma omp for collapse(2) - #endif - //Calculate the standard deviation for each parameter within each considered chain - for (int m = 0; m < Nchains; ++m) - { - for (int j = 0; j < nDraw; ++j) - { - for(int i = 0; i < Ntoys; i++) + if (m == 0) { - StandardDeviation[m][j] += (Draws[m][i][j] - Mean[m][j])*(Draws[m][i][j] - Mean[m][j]); - StandardDeviationFolded[m][j] += (DrawsFolded[m][i][j] - MeanFolded[m][j])*(DrawsFolded[m][i][j] - MeanFolded[m][j]); + StandardDeviationGlobal[j] = 0.0; } - StandardDeviation[m][j] = StandardDeviation[m][j]/(Ntoys-1); - StandardDeviationFolded[m][j] = StandardDeviationFolded[m][j]/(Ntoys-1); - } - } - - #ifdef MULTITHREAD - #pragma omp for - #endif - //Calculate the standard deviation for each parameter combining information from all chains - for (int j = 0; j < nDraw; ++j) - { - for (int m = 0; m < Nchains; ++m) - { StandardDeviationGlobal[j] += StandardDeviation[m][j]; - StandardDeviationGlobalFolded[j] += StandardDeviationFolded[m][j]; } - StandardDeviationGlobal[j] = StandardDeviationGlobal[j]/Nchains; - StandardDeviationGlobalFolded[j] = StandardDeviationGlobalFolded[j]/Nchains; + MeanGlobal[j] = S1_global[j] / (double)TotToys; + StandardDeviationGlobal[j] = StandardDeviationGlobal[j] / (double)Nchains; } #ifdef MULTITHREAD @@ -518,27 +430,24 @@ void CalcRhat() { if(Nchains == 1) { BetweenChainVariance[j] = 0.; - BetweenChainVarianceFolded[j] = 0.; } else { for (int m = 0; m < Nchains; ++m) { - BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]); - BetweenChainVarianceFolded[j] += ( MeanFolded[m][j] - MeanGlobalFolded[j])*( MeanFolded[m][j] - MeanGlobalFolded[j]); + BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]) * Ntoys[m]; } - BetweenChainVariance[j] = BetweenChainVariance[j]*Ntoys/(Nchains-1); - BetweenChainVarianceFolded[j] = BetweenChainVarianceFolded[j]*Ntoys/(Nchains-1); + BetweenChainVariance[j] = BetweenChainVariance[j]/(Nchains-1); } } + int avgNtoys = TotToys/Nchains; #ifdef MULTITHREAD #pragma omp for #endif for (int j = 0; j < nDraw; ++j) { - MarginalPosteriorVariance[j] = (Ntoys-1) * StandardDeviationGlobal[j] /(Ntoys) + BetweenChainVariance[j]/Ntoys; - MarginalPosteriorVarianceFolded[j] = (Ntoys-1) * StandardDeviationGlobalFolded[j] /(Ntoys) + BetweenChainVarianceFolded[j]/Ntoys; + MarginalPosteriorVariance[j] = (avgNtoys-1) * StandardDeviationGlobal[j] / (avgNtoys) + BetweenChainVariance[j]/avgNtoys; } #ifdef MULTITHREAD @@ -548,11 +457,9 @@ void CalcRhat() { for (int j = 0; j < nDraw; ++j) { RHat[j] = sqrt(MarginalPosteriorVariance[j]/StandardDeviationGlobal[j]); - RHatFolded[j] = sqrt(MarginalPosteriorVarianceFolded[j]/StandardDeviationGlobalFolded[j]); //KS: For flat params values can be crazy so cap at 0 CapVariable(RHat[j], 0); - CapVariable(RHatFolded[j], 0); } #ifdef MULTITHREAD @@ -561,12 +468,10 @@ void CalcRhat() { //KS: Additionally calculates effective step size which is an estimate of the sample size required to achieve the same level of precision if that sample was a simple random sample. for (int j = 0; j < nDraw; ++j) { - if(Nchains > 1) EffectiveSampleSize[j] = Nchains * Ntoys * MarginalPosteriorVariance[j] / BetweenChainVariance[j]; - if(Nchains > 1) EffectiveSampleSizeFolded[j] = Nchains * Ntoys * MarginalPosteriorVarianceFolded[j] / BetweenChainVarianceFolded[j]; + if(Nchains > 1) EffectiveSampleSize[j] = TotToys * MarginalPosteriorVariance[j] / BetweenChainVariance[j]; //KS: For flat params values can be crazy so cap at 0 CapVariable(EffectiveSampleSize[j], 0); - CapVariable(EffectiveSampleSizeFolded[j], 0); } #ifdef MULTITHREAD } //End parallel region @@ -612,17 +517,9 @@ void SaveResults() { TH1D *RhatPlot = new TH1D("RhatPlot", "RhatPlot", 200, 0, 2); TH1D *EffectiveSampleSizePlot = new TH1D("EffectiveSampleSizePlot", "EffectiveSampleSizePlot", 400, 0, 10000); - TH1D *StandardDeviationGlobalFoldedPlot = new TH1D("StandardDeviationGlobalFoldedPlot", "StandardDeviationGlobalFoldedPlot", 200, 0, 2); - TH1D *BetweenChainVarianceFoldedPlot = new TH1D("BetweenChainVarianceFoldedPlot", "BetweenChainVarianceFoldedPlot", 200, 0, 2); - TH1D *MarginalPosteriorVarianceFoldedPlot = new TH1D("MarginalPosteriorVarianceFoldedPlot", "MarginalPosteriorVarianceFoldedPlot", 200, 0, 2); - TH1D *RhatFoldedPlot = new TH1D("RhatFoldedPlot", "RhatFoldedPlot", 200, 0, 2); - TH1D *EffectiveSampleSizeFoldedPlot = new TH1D("EffectiveSampleSizeFoldedPlot", "EffectiveSampleSizeFoldedPlot", 400, 0, 10000); - TH1D *RhatLogPlot = new TH1D("RhatLogPlot", "RhatLogPlot", 200, 0, 2); - TH1D *RhatFoldedLogPlot = new TH1D("RhatFoldedLogPlot", "RhatFoldedLogPlot", 200, 0, 2); int Criterium = 0; - int CiteriumFolded = 0; for(int j = 0; j < nDraw; j++) { //KS: Fill only valid parameters @@ -634,26 +531,17 @@ void SaveResults() { RhatPlot->Fill(RHat[j]); EffectiveSampleSizePlot->Fill(EffectiveSampleSize[j]); if(RHat[j] > 1.1) Criterium++; - - - StandardDeviationGlobalFoldedPlot->Fill(StandardDeviationGlobalFolded[j]); - BetweenChainVarianceFoldedPlot->Fill(BetweenChainVarianceFolded[j]); - MarginalPosteriorVarianceFoldedPlot->Fill(MarginalPosteriorVarianceFolded[j]); - RhatFoldedPlot->Fill(RHatFolded[j]); - EffectiveSampleSizeFoldedPlot->Fill(EffectiveSampleSizeFolded[j]); - if(RHatFolded[j] > 1.1) CiteriumFolded++; } else { RhatLogPlot->Fill(RHat[j]); - RhatFoldedLogPlot->Fill(RHatFolded[j]); } } //KS: We set criterium of 1.1 based on Gelman et al. (2003) Bayesian Data Analysis - MACH3LOG_WARN("Number of parameters which has R hat greater than 1.1 is {}({:.2f}%) while for R hat folded {}({:.2f}%)", Criterium, 100*double(Criterium)/double(nDraw), CiteriumFolded, 100*double(CiteriumFolded)/double(nDraw)); + MACH3LOG_WARN("Number of parameters which has R hat greater than 1.1 is {}({:.2f}%)", Criterium, 100*double(Criterium)/double(nDraw)); for(int j = 0; j < nDraw; j++) { - if( (RHat[j] > 1.1 || RHatFolded[j] > 1.1) && ValidPar[j]) + if( (RHat[j] > 1.1) && ValidPar[j]) { MACH3LOG_CRITICAL("Parameter {} has R hat higher than 1.1", BranchNames[j]); } @@ -664,14 +552,7 @@ void SaveResults() { RhatPlot->Write(); EffectiveSampleSizePlot->Write(); - StandardDeviationGlobalFoldedPlot->Write(); - BetweenChainVarianceFoldedPlot->Write(); - MarginalPosteriorVarianceFoldedPlot->Write(); - RhatFoldedPlot->Write(); - EffectiveSampleSizeFoldedPlot->Write(); - RhatLogPlot->Write(); - RhatFoldedLogPlot->Write(); //KS: Now we make fancy canvases, consider some function to have less copy pasting auto TempCanvas = std::make_unique("Canvas", "Canvas", 1024, 1024); @@ -686,8 +567,6 @@ void SaveResults() { RhatPlot->GetXaxis()->SetTitle("R hat"); RhatPlot->SetLineColor(kRed); RhatPlot->SetFillColor(kRed); - RhatFoldedPlot->SetLineColor(kBlue); - RhatFoldedPlot->SetFillColor(kBlue); TLegend *Legend = new TLegend(0.55, 0.6, 0.9, 0.9); Legend->SetTextSize(0.04); @@ -698,10 +577,8 @@ void SaveResults() { Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); Legend->AddEntry(RhatPlot, "Rhat Gelman 2013", "l"); - Legend->AddEntry(RhatFoldedPlot, "Rhat-Folded Gelman 2021", "l"); RhatPlot->Draw(); - RhatFoldedPlot->Draw("same"); Legend->Draw("same"); TempCanvas->Write("Rhat"); delete Legend; @@ -711,8 +588,6 @@ void SaveResults() { RhatLogPlot->GetXaxis()->SetTitle("R hat for LogL"); RhatLogPlot->SetLineColor(kRed); RhatLogPlot->SetFillColor(kRed); - RhatFoldedLogPlot->SetLineColor(kBlue); - RhatFoldedLogPlot->SetFillColor(kBlue); Legend = new TLegend(0.55, 0.6, 0.9, 0.9); Legend->SetTextSize(0.04); @@ -723,10 +598,8 @@ void SaveResults() { Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); Legend->AddEntry(RhatLogPlot, "Rhat Gelman 2013", "l"); - Legend->AddEntry(RhatFoldedLogPlot, "Rhat-Folded Gelman 2021", "l"); RhatLogPlot->Draw(); - RhatFoldedLogPlot->Draw("same"); Legend->Draw("same"); TempCanvas->Write("RhatLog"); delete Legend; @@ -735,7 +608,6 @@ void SaveResults() { //Now canvas for effective sample size EffectiveSampleSizePlot->GetXaxis()->SetTitle("S_{eff, BDA2}"); EffectiveSampleSizePlot->SetLineColor(kRed); - EffectiveSampleSizeFoldedPlot->SetLineColor(kBlue); Legend = new TLegend(0.45, 0.6, 0.9, 0.9); Legend->SetTextSize(0.03); @@ -746,15 +618,11 @@ void SaveResults() { const double Mean1 = EffectiveSampleSizePlot->GetMean(); const double RMS1 = EffectiveSampleSizePlot->GetRMS(); - const double Mean2 = EffectiveSampleSizeFoldedPlot->GetMean(); - const double RMS2 = EffectiveSampleSizeFoldedPlot->GetRMS(); Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); Legend->AddEntry(EffectiveSampleSizePlot, Form("S_{eff, BDA2} #mu = %.2f, #sigma = %.2f",Mean1 ,RMS1), "l"); - Legend->AddEntry(EffectiveSampleSizeFoldedPlot, Form("S_{eff, BDA2} Folded, #mu = %.2f, #sigma = %.2f",Mean2 ,RMS2), "l"); EffectiveSampleSizePlot->Draw(); - EffectiveSampleSizeFoldedPlot->Draw("same"); Legend->Draw("same"); TempCanvas->Write("EffectiveSampleSize"); @@ -765,16 +633,9 @@ void SaveResults() { delete RhatPlot; delete EffectiveSampleSizePlot; - delete StandardDeviationGlobalFoldedPlot; - delete BetweenChainVarianceFoldedPlot; - delete MarginalPosteriorVarianceFoldedPlot; - delete RhatFoldedPlot; - delete EffectiveSampleSizeFoldedPlot; - delete Legend; delete RhatLogPlot; - delete RhatFoldedLogPlot; DiagFile->Close(); delete DiagFile; @@ -795,36 +656,21 @@ void DestroyArrays() { delete[] RHat; delete[] EffectiveSampleSize; - delete[] MeanGlobalFolded; - delete[] StandardDeviationGlobalFolded; - delete[] BetweenChainVarianceFolded; - delete[] MarginalPosteriorVarianceFolded; - delete[] RHatFolded; - delete[] EffectiveSampleSizeFolded; - for(int m = 0; m < Nchains; m++) { - for(int i = 0; i < Ntoys; i++) - { - delete[] Draws[m][i]; - delete[] DrawsFolded[m][i]; - } - delete[] Draws[m]; delete[] Mean[m]; delete[] StandardDeviation[m]; - - delete[] DrawsFolded[m]; - delete[] MeanFolded[m]; - delete[] StandardDeviationFolded[m]; + delete[] S1_chain[m]; + delete[] S2_chain[m]; } - delete[] Draws; delete[] Mean; delete[] StandardDeviation; + delete[] S1_chain; + delete[] S2_chain; + delete[] S1_global; + delete[] S2_global; - delete[] DrawsFolded; - delete[] MedianArr; - delete[] MeanFolded; - delete[] StandardDeviationFolded; + delete[] Ntoys; } // ******************* diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp new file mode 100644 index 000000000..a0f9976d2 --- /dev/null +++ b/Diagnostics/RHat_HighMem.cpp @@ -0,0 +1,845 @@ +// ROOT includes +#include "TObjArray.h" +#include "TChain.h" +#include "TFile.h" +#include "TBranch.h" +#include "TCanvas.h" +#include "TLine.h" +#include "TLegend.h" +#include "TString.h" +#include "TH1.h" +#include "TRandom3.h" +#include "TStopwatch.h" +#include "TColor.h" +#include "TStyle.h" +#include "TROOT.h" + +#ifdef MULTITHREAD +#include "omp.h" +#endif + +// MaCh3 includes +#include "manager/manager.h" + +/// @file RHat.cpp +/// @brief This executable calculates the \f$ \hat{R} \f$ estimator for Markov Chain Monte Carlo (MCMC) convergence. +/// +/// KS: This exe is meant to calculate the \f$ \hat{R} \f$ estimator. For a well-converged chain, this distribution +/// should be centered at one. The \f$ \hat{R} \f$ statistic is used to assess the convergence of MCMC simulations +/// and helps determine whether the chains have reached a stable distribution. +/// +/// @cite gelman2019. + +// ******************* +int Ntoys; +int Nchains; + +int nDraw; + +std::vector BranchNames; +std::vector MCMCFile; +std::vector ValidPar; + +double ***Draws; + +double** Mean; +double** StandardDeviation; + +double* MeanGlobal; +double* StandardDeviationGlobal; + +double* BetweenChainVariance; +double* MarginalPosteriorVariance; +double* RHat; +double* EffectiveSampleSize; + +double ***DrawsFolded; +double* MedianArr; + +double** MeanFolded; +double** StandardDeviationFolded; + +double* MeanGlobalFolded; +double* StandardDeviationGlobalFolded; + +double* BetweenChainVarianceFolded; +double* MarginalPosteriorVarianceFolded; +double* RHatFolded; +double* EffectiveSampleSizeFolded; +// ******************* +void PrepareChains(); +void InitialiseArrays(); + +void RunDiagnostic(); +void CalcRhat(); + +void SaveResults(); +void DestroyArrays(); +double CalcMedian(double arr[], int size); + +void CapVariable(double var, double cap); + +// ******************* +int main(int argc, char *argv[]) { +// ******************* + + SetMaCh3LoggerFormat(); + MaCh3Utils::MaCh3Welcome(); + + Draws = nullptr; + Mean = nullptr; + StandardDeviation = nullptr; + + MeanGlobal = nullptr; + StandardDeviationGlobal = nullptr; + + BetweenChainVariance = nullptr; + MarginalPosteriorVariance = nullptr; + RHat = nullptr; + EffectiveSampleSize = nullptr; + + DrawsFolded = nullptr; + MedianArr = nullptr; + MeanFolded = nullptr; + StandardDeviationFolded = nullptr; + + MeanGlobalFolded = nullptr; + StandardDeviationGlobalFolded = nullptr; + + BetweenChainVarianceFolded = nullptr; + MarginalPosteriorVarianceFolded = nullptr; + RHatFolded = nullptr; + EffectiveSampleSizeFolded = nullptr; + + Nchains = 0; + + if (argc == 1 || argc == 2) + { + MACH3LOG_ERROR("Wrong arguments"); + MACH3LOG_ERROR("./RHat Ntoys MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + Ntoys = atoi(argv[1]); + + //KS Gelman suggests to diagnose on more than one chain + for (int i = 2; i < argc; i++) + { + MCMCFile.push_back(std::string(argv[i])); + MACH3LOG_INFO("Adding file: {}", MCMCFile.back()); + Nchains++; + } + + if(Ntoys < 1) + { + MACH3LOG_ERROR("You specified {} specify larger greater than 0", Ntoys); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + if(Nchains == 1) + { + MACH3LOG_WARN("Gelman is going to be sad :(. He suggested you should use more than one chain (at least 4). Code works fine for one chain, however, estimator might be biased."); + MACH3LOG_WARN("Multiple chains are more likely to reveal multimodality and poor adaptation or mixing:"); + } + MACH3LOG_INFO("Diagnosing {} chains, with {} toys", Nchains, Ntoys); + + PrepareChains(); + + InitialiseArrays(); + + //KS: Main function + RunDiagnostic(); + + SaveResults(); + + DestroyArrays(); + + return 0; +} + +// ******************* +// Load chain and prepare toys +void PrepareChains() { +// ******************* + auto rnd = std::make_unique(0); + + MACH3LOG_INFO("Generating {}", Ntoys); + + TStopwatch clock; + clock.Start(); + + std::vector BurnIn(Nchains); + std::vector nEntries(Nchains); + std::vector nBranches(Nchains); + std::vector step(Nchains); + + Draws = new double**[Nchains](); + DrawsFolded = new double**[Nchains](); + + // KS: This can reduce time necessary for caching even by half + #ifdef MULTITHREAD + //ROOT::EnableImplicitMT(); + #endif + + // Open the Chain + //It is tempting to multithread here but unfortunately, ROOT files are not thread safe :( + for (int m = 0; m < Nchains; m++) + { + TChain* Chain = new TChain("posteriors"); + Chain->Add(MCMCFile[m].c_str()); + MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); + nEntries[m] = int(Chain->GetEntries()); + + // Set the step cut to be 20% + BurnIn[m] = nEntries[m]/5; + + // Get the list of branches + TObjArray* brlis = Chain->GetListOfBranches(); + + // Get the number of branches + nBranches[m] = brlis->GetEntries(); + + if(m == 0) BranchNames.reserve(nBranches[m]); + + // Set all the branches to off + Chain->SetBranchStatus("*", false); + + // Loop over the number of branches + // Find the name and how many of each systematic we have + for (int i = 0; i < nBranches[m]; i++) + { + // Get the TBranch and its name + TBranch* br = static_cast(brlis->At(i)); + if(!br){ + MACH3LOG_ERROR("Invalid branch at position {}", i); + throw MaCh3Exception(__FILE__,__LINE__); + } + TString bname = br->GetName(); + + // Read in the step + if (bname == "step") { + Chain->SetBranchStatus(bname, true); + Chain->SetBranchAddress(bname, &step[m]); + } + //Count all branches + else if (bname.BeginsWith("PCA_") || bname.BeginsWith("accProb") || bname.BeginsWith("stepTime") ) + { + continue; + } + else + { + //KS: Save branch name only for one chain, we assume all chains have the same branches, otherwise this doesn't make sense either way + if(m == 0) + { + BranchNames.push_back(bname); + //KS: We calculate R Hat also for LogL, just in case, however we plot them separately + if(bname.BeginsWith("LogL")) + { + ValidPar.push_back(false); + } + else + { + ValidPar.push_back(true); + } + } + Chain->SetBranchStatus(bname, true); + MACH3LOG_DEBUG("{}", bname); + } + } + + if(m == 0) nDraw = int(BranchNames.size()); + + //TN: Qualitatively faster sanity check, with the very same outcome (all chains have the same #branches) + if(m > 0) + { + if(nBranches[m] != nBranches[0]) + { + MACH3LOG_ERROR("Ups, something went wrong, chain {} called {} has {} branches, while 0 called {} has {} branches", m, MCMCFile[m], nBranches[m], MCMCFile[0], nBranches[0]); + MACH3LOG_ERROR("All chains should have the same number of branches"); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + } + + //TN: move the Draws here, so we need to iterate over every chain only once + Draws[m] = new double*[Ntoys](); + DrawsFolded[m] = new double*[Ntoys](); + for(int i = 0; i < Ntoys; i++) + { + Draws[m][i] = new double[nDraw](); + DrawsFolded[m][i] = new double[nDraw](); + for(int j = 0; j < nDraw; j++) + { + Draws[m][i][j] = 0.; + DrawsFolded[m][i][j] = 0.; + } + } + + //TN: move looping over toys here, so we don't need to loop over chains more than once + if(BurnIn[m] >= nEntries[m]) + { + MACH3LOG_ERROR("You are running on a chain shorter than BurnIn cut"); + MACH3LOG_ERROR("Number of entries {} BurnIn cut {}", nEntries[m], BurnIn[m]); + MACH3LOG_ERROR("You will run into the infinite loop"); + MACH3LOG_ERROR("You can make a new chain or modify BurnIn cut"); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + for (int i = 0; i < Ntoys; i++) + { + // Get a random entry after burn in + int entry = int(nEntries[m]*rnd->Rndm()); + + Chain->GetEntry(entry); + + // If we have combined chains by hadd need to check the step in the chain + // Note, entry is not necessarily the same as the step due to merged ROOT files, so can't choose an entry in the range BurnIn - nEntries :( + if (step[m] < BurnIn[m]) + { + i--; + continue; + } + + // Output some info for the user + if (Ntoys > 10 && i % (Ntoys/10) == 0) { + MaCh3Utils::PrintProgressBar(i+m*Ntoys, static_cast(Ntoys)*Nchains); + MACH3LOG_DEBUG("Getting random entry {}", entry); + } + + // Set the branch addresses for params + for (int j = 0; j < nDraw; ++j) + { + Chain->SetBranchAddress(BranchNames[j].Data(), &Draws[m][i][j]); + } + Chain->GetEntry(entry); + + }//end loop over toys + + //TN: There, we now don't need to keep the chain in memory anymore + delete Chain; + } + + //KS: Now prepare folded draws, quoting Gelman + //"We propose to report the maximum of rank normalized split-Rb and rank normalized folded-split-Rb for each parameter" + MedianArr = new double[nDraw](); + #ifdef MULTITHREAD + #pragma omp parallel for + #endif + for(int j = 0; j < nDraw; j++) + { + MedianArr[j] = 0.; + std::vector TempDraws(static_cast(Ntoys) * Nchains); + for(int m = 0; m < Nchains; m++) + { + for(int i = 0; i < Ntoys; i++) + { + const int im = i+m; + TempDraws[im] = Draws[m][i][j]; + } + } + MedianArr[j] = CalcMedian(TempDraws.data(), Ntoys*Nchains); + } + + #ifdef MULTITHREAD + #pragma omp parallel for collapse(3) + #endif + for(int m = 0; m < Nchains; m++) + { + for(int i = 0; i < Ntoys; i++) + { + for(int j = 0; j < nDraw; j++) + { + DrawsFolded[m][i][j] = std::fabs(Draws[m][i][j] - MedianArr[j]); + } + } + } + clock.Stop(); + MACH3LOG_INFO("Finished calculating Toys, it took {:.2f}s to finish", clock.RealTime()); +} + +// ******************* +// Create all arrays we are going to use later +void InitialiseArrays() { +// ******************* + + MACH3LOG_INFO("Initialising arrays"); + Mean = new double*[Nchains](); + StandardDeviation = new double*[Nchains](); + + MeanGlobal = new double[nDraw](); + StandardDeviationGlobal = new double[nDraw](); + BetweenChainVariance = new double[nDraw](); + + MarginalPosteriorVariance = new double[nDraw](); + RHat = new double[nDraw](); + EffectiveSampleSize = new double[nDraw](); + + MeanFolded = new double*[Nchains](); + StandardDeviationFolded = new double*[Nchains](); + + MeanGlobalFolded = new double[nDraw](); + StandardDeviationGlobalFolded = new double[nDraw](); + BetweenChainVarianceFolded = new double[nDraw](); + + MarginalPosteriorVarianceFolded = new double[nDraw](); + RHatFolded = new double[nDraw](); + EffectiveSampleSizeFolded = new double[nDraw](); + + for (int m = 0; m < Nchains; ++m) + { + Mean[m] = new double[nDraw](); + StandardDeviation[m] = new double[nDraw](); + + MeanFolded[m] = new double[nDraw](); + StandardDeviationFolded[m] = new double[nDraw](); + for (int j = 0; j < nDraw; ++j) + { + Mean[m][j] = 0.; + StandardDeviation[m][j] = 0.; + + MeanFolded[m][j] = 0.; + StandardDeviationFolded[m][j] = 0.; + if(m == 0) + { + MeanGlobal[j] = 0.; + StandardDeviationGlobal[j] = 0.; + BetweenChainVariance[j] = 0.; + MarginalPosteriorVariance[j] = 0.; + RHat[j] = 0.; + EffectiveSampleSize[j] = 0.; + + MeanGlobalFolded[j] = 0.; + StandardDeviationGlobalFolded[j] = 0.; + BetweenChainVarianceFolded[j] = 0.; + MarginalPosteriorVarianceFolded[j] = 0.; + RHatFolded[j] = 0.; + EffectiveSampleSizeFolded[j] = 0.; + } + } + } +} + +// ******************* +void RunDiagnostic() { +// ******************* + CalcRhat(); + //In case in future we expand this +} + +// ******************* +//KS: Based on Gelman et. al. arXiv:1903.08008v5 +// Probably most of it could be moved cleverly to MCMC Processor, keep it separate for now +void CalcRhat() { +// ******************* + + TStopwatch clock; + clock.Start(); + + //KS: Start parallel region + // If we would like to do this for thousands of chains we might consider using GPU for this + #ifdef MULTITHREAD + #pragma omp parallel + { + #endif + + #ifdef MULTITHREAD + #pragma omp for collapse(2) + #endif + //KS: loop over chains and draws are independent so might as well collapse for sweet cache hits + //Calculate the mean for each parameter within each considered chain + for (int m = 0; m < Nchains; ++m) + { + for (int j = 0; j < nDraw; ++j) + { + for(int i = 0; i < Ntoys; i++) + { + Mean[m][j] += Draws[m][i][j]; + MeanFolded[m][j] += DrawsFolded[m][i][j]; + } + Mean[m][j] = Mean[m][j]/Ntoys; + MeanFolded[m][j] = MeanFolded[m][j]/Ntoys; + } + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + //Calculate the mean for each parameter global means we include information from several chains + for (int j = 0; j < nDraw; ++j) + { + for (int m = 0; m < Nchains; ++m) + { + MeanGlobal[j] += Mean[m][j]; + MeanGlobalFolded[j] += MeanFolded[m][j]; + } + MeanGlobal[j] = MeanGlobal[j]/Nchains; + MeanGlobalFolded[j] = MeanGlobalFolded[j]/Nchains; + } + + + #ifdef MULTITHREAD + #pragma omp for collapse(2) + #endif + //Calculate the standard deviation for each parameter within each considered chain + for (int m = 0; m < Nchains; ++m) + { + for (int j = 0; j < nDraw; ++j) + { + for(int i = 0; i < Ntoys; i++) + { + StandardDeviation[m][j] += (Draws[m][i][j] - Mean[m][j])*(Draws[m][i][j] - Mean[m][j]); + StandardDeviationFolded[m][j] += (DrawsFolded[m][i][j] - MeanFolded[m][j])*(DrawsFolded[m][i][j] - MeanFolded[m][j]); + } + StandardDeviation[m][j] = StandardDeviation[m][j]/(Ntoys-1); + StandardDeviationFolded[m][j] = StandardDeviationFolded[m][j]/(Ntoys-1); + } + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + //Calculate the standard deviation for each parameter combining information from all chains + for (int j = 0; j < nDraw; ++j) + { + for (int m = 0; m < Nchains; ++m) + { + StandardDeviationGlobal[j] += StandardDeviation[m][j]; + StandardDeviationGlobalFolded[j] += StandardDeviationFolded[m][j]; + } + StandardDeviationGlobal[j] = StandardDeviationGlobal[j]/Nchains; + StandardDeviationGlobalFolded[j] = StandardDeviationGlobalFolded[j]/Nchains; + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + for (int j = 0; j < nDraw; ++j) + { + //KS: This term only makes sense if we have at least 2 chains + if(Nchains == 1) + { + BetweenChainVariance[j] = 0.; + BetweenChainVarianceFolded[j] = 0.; + } + else + { + for (int m = 0; m < Nchains; ++m) + { + BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]); + BetweenChainVarianceFolded[j] += ( MeanFolded[m][j] - MeanGlobalFolded[j])*( MeanFolded[m][j] - MeanGlobalFolded[j]); + } + BetweenChainVariance[j] = BetweenChainVariance[j]*Ntoys/(Nchains-1); + BetweenChainVarianceFolded[j] = BetweenChainVarianceFolded[j]*Ntoys/(Nchains-1); + } + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + for (int j = 0; j < nDraw; ++j) + { + MarginalPosteriorVariance[j] = (Ntoys-1) * StandardDeviationGlobal[j] /(Ntoys) + BetweenChainVariance[j]/Ntoys; + MarginalPosteriorVarianceFolded[j] = (Ntoys-1) * StandardDeviationGlobalFolded[j] /(Ntoys) + BetweenChainVarianceFolded[j]/Ntoys; + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + //Finally calculate our estimator + for (int j = 0; j < nDraw; ++j) + { + RHat[j] = sqrt(MarginalPosteriorVariance[j]/StandardDeviationGlobal[j]); + RHatFolded[j] = sqrt(MarginalPosteriorVarianceFolded[j]/StandardDeviationGlobalFolded[j]); + + //KS: For flat params values can be crazy so cap at 0 + CapVariable(RHat[j], 0); + CapVariable(RHatFolded[j], 0); + } + + #ifdef MULTITHREAD + #pragma omp for + #endif + //KS: Additionally calculates effective step size which is an estimate of the sample size required to achieve the same level of precision if that sample was a simple random sample. + for (int j = 0; j < nDraw; ++j) + { + if(Nchains > 1) EffectiveSampleSize[j] = Nchains * Ntoys * MarginalPosteriorVariance[j] / BetweenChainVariance[j]; + if(Nchains > 1) EffectiveSampleSizeFolded[j] = Nchains * Ntoys * MarginalPosteriorVarianceFolded[j] / BetweenChainVarianceFolded[j]; + + //KS: For flat params values can be crazy so cap at 0 + CapVariable(EffectiveSampleSize[j], 0); + CapVariable(EffectiveSampleSizeFolded[j], 0); + } + #ifdef MULTITHREAD + } //End parallel region + #endif + + clock.Stop(); + MACH3LOG_INFO("Finished calculating RHat, it took {:.2f}s to finish", clock.RealTime()); +} + + +// ******************* +void SaveResults() { +// ******************* + #pragma GCC diagnostic ignored "-Wfloat-conversion" + + std::string NameTemp = ""; + //KS: If we run over many many chains there is danger that name will be so absurdly long we run over system limit and job will be killed :( + if(Nchains < 5) + { + for (int i = 0; i < Nchains; i++) + { + std::string temp = MCMCFile[i]; + + while (temp.find(".root") != std::string::npos) { + temp = temp.substr(0, temp.find(".root")); + } + + NameTemp = NameTemp + temp + "_"; + } + } + else { + NameTemp = std::to_string(Nchains) + "Chains" + "_"; + } + NameTemp += "diag.root"; + + TFile* DiagFile = new TFile(NameTemp.c_str(), "recreate"); + + DiagFile->cd(); + + TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", 200, 0, 2); + TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", 200, 0, 2); + TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", 200, 0, 2); + TH1D *RhatPlot = new TH1D("RhatPlot", "RhatPlot", 200, 0, 2); + TH1D *EffectiveSampleSizePlot = new TH1D("EffectiveSampleSizePlot", "EffectiveSampleSizePlot", 400, 0, 10000); + + TH1D *StandardDeviationGlobalFoldedPlot = new TH1D("StandardDeviationGlobalFoldedPlot", "StandardDeviationGlobalFoldedPlot", 200, 0, 2); + TH1D *BetweenChainVarianceFoldedPlot = new TH1D("BetweenChainVarianceFoldedPlot", "BetweenChainVarianceFoldedPlot", 200, 0, 2); + TH1D *MarginalPosteriorVarianceFoldedPlot = new TH1D("MarginalPosteriorVarianceFoldedPlot", "MarginalPosteriorVarianceFoldedPlot", 200, 0, 2); + TH1D *RhatFoldedPlot = new TH1D("RhatFoldedPlot", "RhatFoldedPlot", 200, 0, 2); + TH1D *EffectiveSampleSizeFoldedPlot = new TH1D("EffectiveSampleSizeFoldedPlot", "EffectiveSampleSizeFoldedPlot", 400, 0, 10000); + + TH1D *RhatLogPlot = new TH1D("RhatLogPlot", "RhatLogPlot", 200, 0, 2); + TH1D *RhatFoldedLogPlot = new TH1D("RhatFoldedLogPlot", "RhatFoldedLogPlot", 200, 0, 2); + + int Criterium = 0; + int CiteriumFolded = 0; + for(int j = 0; j < nDraw; j++) + { + //KS: Fill only valid parameters + if(ValidPar[j]) + { + StandardDeviationGlobalPlot->Fill(StandardDeviationGlobal[j]); + BetweenChainVariancePlot->Fill(BetweenChainVariance[j]); + MarginalPosteriorVariancePlot->Fill(MarginalPosteriorVariance[j]); + RhatPlot->Fill(RHat[j]); + EffectiveSampleSizePlot->Fill(EffectiveSampleSize[j]); + if(RHat[j] > 1.1) Criterium++; + + + StandardDeviationGlobalFoldedPlot->Fill(StandardDeviationGlobalFolded[j]); + BetweenChainVarianceFoldedPlot->Fill(BetweenChainVarianceFolded[j]); + MarginalPosteriorVarianceFoldedPlot->Fill(MarginalPosteriorVarianceFolded[j]); + RhatFoldedPlot->Fill(RHatFolded[j]); + EffectiveSampleSizeFoldedPlot->Fill(EffectiveSampleSizeFolded[j]); + if(RHatFolded[j] > 1.1) CiteriumFolded++; + } + else + { + RhatLogPlot->Fill(RHat[j]); + RhatFoldedLogPlot->Fill(RHatFolded[j]); + } + } + //KS: We set criterium of 1.1 based on Gelman et al. (2003) Bayesian Data Analysis + MACH3LOG_WARN("Number of parameters which has R hat greater than 1.1 is {}({:.2f}%) while for R hat folded {}({:.2f}%)", Criterium, 100*double(Criterium)/double(nDraw), CiteriumFolded, 100*double(CiteriumFolded)/double(nDraw)); + for(int j = 0; j < nDraw; j++) + { + if( (RHat[j] > 1.1 || RHatFolded[j] > 1.1) && ValidPar[j]) + { + MACH3LOG_CRITICAL("Parameter {} has R hat higher than 1.1", BranchNames[j]); + } + } + StandardDeviationGlobalPlot->Write(); + BetweenChainVariancePlot->Write(); + MarginalPosteriorVariancePlot->Write(); + RhatPlot->Write(); + EffectiveSampleSizePlot->Write(); + + StandardDeviationGlobalFoldedPlot->Write(); + BetweenChainVarianceFoldedPlot->Write(); + MarginalPosteriorVarianceFoldedPlot->Write(); + RhatFoldedPlot->Write(); + EffectiveSampleSizeFoldedPlot->Write(); + + RhatLogPlot->Write(); + RhatFoldedLogPlot->Write(); + + //KS: Now we make fancy canvases, consider some function to have less copy pasting + auto TempCanvas = std::make_unique("Canvas", "Canvas", 1024, 1024); + gStyle->SetOptStat(0); + TempCanvas->SetGridx(); + TempCanvas->SetGridy(); + + // Random line to write useful information to TLegend + auto TempLine = std::make_unique(0, 0, 0, 0); + TempLine->SetLineColor(kBlack); + + RhatPlot->GetXaxis()->SetTitle("R hat"); + RhatPlot->SetLineColor(kRed); + RhatPlot->SetFillColor(kRed); + RhatFoldedPlot->SetLineColor(kBlue); + RhatFoldedPlot->SetFillColor(kBlue); + + TLegend *Legend = new TLegend(0.55, 0.6, 0.9, 0.9); + Legend->SetTextSize(0.04); + Legend->SetFillColor(0); + Legend->SetFillStyle(0); + Legend->SetLineWidth(0); + Legend->SetLineColor(0); + + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(RhatPlot, "Rhat Gelman 2013", "l"); + Legend->AddEntry(RhatFoldedPlot, "Rhat-Folded Gelman 2021", "l"); + + RhatPlot->Draw(); + RhatFoldedPlot->Draw("same"); + Legend->Draw("same"); + TempCanvas->Write("Rhat"); + delete Legend; + Legend = nullptr; + + //Now R hat for log L + RhatLogPlot->GetXaxis()->SetTitle("R hat for LogL"); + RhatLogPlot->SetLineColor(kRed); + RhatLogPlot->SetFillColor(kRed); + RhatFoldedLogPlot->SetLineColor(kBlue); + RhatFoldedLogPlot->SetFillColor(kBlue); + + Legend = new TLegend(0.55, 0.6, 0.9, 0.9); + Legend->SetTextSize(0.04); + Legend->SetFillColor(0); + Legend->SetFillStyle(0); + Legend->SetLineWidth(0); + Legend->SetLineColor(0); + + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(RhatLogPlot, "Rhat Gelman 2013", "l"); + Legend->AddEntry(RhatFoldedLogPlot, "Rhat-Folded Gelman 2021", "l"); + + RhatLogPlot->Draw(); + RhatFoldedLogPlot->Draw("same"); + Legend->Draw("same"); + TempCanvas->Write("RhatLog"); + delete Legend; + Legend = nullptr; + + //Now canvas for effective sample size + EffectiveSampleSizePlot->GetXaxis()->SetTitle("S_{eff, BDA2}"); + EffectiveSampleSizePlot->SetLineColor(kRed); + EffectiveSampleSizeFoldedPlot->SetLineColor(kBlue); + + Legend = new TLegend(0.45, 0.6, 0.9, 0.9); + Legend->SetTextSize(0.03); + Legend->SetFillColor(0); + Legend->SetFillStyle(0); + Legend->SetLineWidth(0); + Legend->SetLineColor(0); + + const double Mean1 = EffectiveSampleSizePlot->GetMean(); + const double RMS1 = EffectiveSampleSizePlot->GetRMS(); + const double Mean2 = EffectiveSampleSizeFoldedPlot->GetMean(); + const double RMS2 = EffectiveSampleSizeFoldedPlot->GetRMS(); + + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(EffectiveSampleSizePlot, Form("S_{eff, BDA2} #mu = %.2f, #sigma = %.2f",Mean1 ,RMS1), "l"); + Legend->AddEntry(EffectiveSampleSizeFoldedPlot, Form("S_{eff, BDA2} Folded, #mu = %.2f, #sigma = %.2f",Mean2 ,RMS2), "l"); + + EffectiveSampleSizePlot->Draw(); + EffectiveSampleSizeFoldedPlot->Draw("same"); + Legend->Draw("same"); + TempCanvas->Write("EffectiveSampleSize"); + + //Fancy memory cleaning + delete StandardDeviationGlobalPlot; + delete BetweenChainVariancePlot; + delete MarginalPosteriorVariancePlot; + delete RhatPlot; + delete EffectiveSampleSizePlot; + + delete StandardDeviationGlobalFoldedPlot; + delete BetweenChainVarianceFoldedPlot; + delete MarginalPosteriorVarianceFoldedPlot; + delete RhatFoldedPlot; + delete EffectiveSampleSizeFoldedPlot; + + delete Legend; + + delete RhatLogPlot; + delete RhatFoldedLogPlot; + + DiagFile->Close(); + delete DiagFile; + + MACH3LOG_INFO("Finished and wrote results to {}", NameTemp); +} + +// ******************* +//KS: Pseudo destructor +void DestroyArrays() { +// ******************* + + MACH3LOG_INFO("Killing all arrays"); + delete[] MeanGlobal; + delete[] StandardDeviationGlobal; + delete[] BetweenChainVariance; + delete[] MarginalPosteriorVariance; + delete[] RHat; + delete[] EffectiveSampleSize; + + delete[] MeanGlobalFolded; + delete[] StandardDeviationGlobalFolded; + delete[] BetweenChainVarianceFolded; + delete[] MarginalPosteriorVarianceFolded; + delete[] RHatFolded; + delete[] EffectiveSampleSizeFolded; + + for(int m = 0; m < Nchains; m++) + { + for(int i = 0; i < Ntoys; i++) + { + delete[] Draws[m][i]; + delete[] DrawsFolded[m][i]; + } + delete[] Draws[m]; + delete[] Mean[m]; + delete[] StandardDeviation[m]; + + delete[] DrawsFolded[m]; + delete[] MeanFolded[m]; + delete[] StandardDeviationFolded[m]; + } + delete[] Draws; + delete[] Mean; + delete[] StandardDeviation; + + delete[] DrawsFolded; + delete[] MedianArr; + delete[] MeanFolded; + delete[] StandardDeviationFolded; +} + +// ******************* +//calculate median +double CalcMedian(double arr[], const int size) { +// ******************* + std::sort(arr, arr+size); + if (size % 2 != 0) + return arr[size/2]; + return (arr[(size-1)/2] + arr[size/2])/2.0; +} + +// ******************* +//calculate median +void CapVariable(double var, const double cap) { +// ******************* + if(std::isnan(var) || !std::isfinite(var)) var = cap; +} From 21add2d61731f2f32d5f89dbe99ba4196535e05b Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 14:18:46 -0500 Subject: [PATCH 02/13] Added new RHat calculator to CMakeLists --- Diagnostics/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Diagnostics/CMakeLists.txt b/Diagnostics/CMakeLists.txt index 63e0ee838..cae39a0ff 100644 --- a/Diagnostics/CMakeLists.txt +++ b/Diagnostics/CMakeLists.txt @@ -3,6 +3,7 @@ add_custom_target(DiagApps) foreach(app DiagMCMC RHat + RHat_HighMem GetPenaltyTerm ProcessMCMC CombineMaCh3Chains From 3278fa1a0acc6b5d2667e25afe6ef2e7c6cce808 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 14:19:16 -0500 Subject: [PATCH 03/13] Fixed several compiler errors for new RHat calculator --- Diagnostics/RHat.cpp | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 3be941df5..7d9776633 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -106,8 +106,6 @@ int main(int argc, char *argv[]) { throw MaCh3Exception(__FILE__ , __LINE__ ); } - Ntoys = atoi(argv[1]); - //KS Gelman suggests to diagnose on more than one chain for (int i = 1; i < argc; i++) { @@ -172,7 +170,7 @@ void PrepareChains() { TotToys += Ntoys[m]; MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); - MACH3LOG_INFO("Generating {}", Ntoys); + MACH3LOG_INFO("Generating {}", Ntoys[m]); // Set the step cut to be 20% BurnIn[m] = nEntries[m]/5; @@ -298,8 +296,8 @@ void PrepareChains() { } // Output some info for the user - if (Ntoys > 10 && i % (Ntoys/10) == 0) { - MaCh3Utils::PrintProgressBar(i+m*Ntoys, static_cast(Ntoys)*Nchains); + if (Ntoys[m] > 10 && i % (Ntoys[m]/10) == 0) { + MaCh3Utils::PrintProgressBar(i+m*Ntoys[m], static_cast(Ntoys[m])*Nchains); MACH3LOG_DEBUG("Getting random entry {}", entry); } @@ -398,8 +396,8 @@ void CalcRhat() { { for (int j = 0; j < nDraw; ++j) { - Mean[m][j] += S1_chain[m][j] / (double)Ntoys[m]; - StandardDeviation[m][j] = S2_chain[m][j]/(double)Ntoys[m] - Mean[m][j]*Mean[m][j]; + Mean[m][j] += S1_chain[m][j] / static_cast(Ntoys[m]); + StandardDeviation[m][j] = S2_chain[m][j]/static_cast(Ntoys[m]) - Mean[m][j]*Mean[m][j]; } } @@ -417,8 +415,8 @@ void CalcRhat() { } StandardDeviationGlobal[j] += StandardDeviation[m][j]; } - MeanGlobal[j] = S1_global[j] / (double)TotToys; - StandardDeviationGlobal[j] = StandardDeviationGlobal[j] / (double)Nchains; + MeanGlobal[j] = S1_global[j] / static_cast(TotToys); + StandardDeviationGlobal[j] = StandardDeviationGlobal[j] / static_cast(Nchains); } #ifdef MULTITHREAD @@ -575,7 +573,7 @@ void SaveResults() { Legend->SetLineWidth(0); Legend->SetLineColor(0); - Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), ""); Legend->AddEntry(RhatPlot, "Rhat Gelman 2013", "l"); RhatPlot->Draw(); @@ -596,7 +594,7 @@ void SaveResults() { Legend->SetLineWidth(0); Legend->SetLineColor(0); - Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), ""); Legend->AddEntry(RhatLogPlot, "Rhat Gelman 2013", "l"); RhatLogPlot->Draw(); @@ -619,7 +617,7 @@ void SaveResults() { const double Mean1 = EffectiveSampleSizePlot->GetMean(); const double RMS1 = EffectiveSampleSizePlot->GetRMS(); - Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", Ntoys, Nchains), ""); + Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), ""); Legend->AddEntry(EffectiveSampleSizePlot, Form("S_{eff, BDA2} #mu = %.2f, #sigma = %.2f",Mean1 ,RMS1), "l"); EffectiveSampleSizePlot->Draw(); From 28eda9ee33bff86fd38e246e723813234e6744db Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 14:23:56 -0500 Subject: [PATCH 04/13] Small comment update --- Diagnostics/RHat.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 7d9776633..1e6d5398f 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -392,6 +392,7 @@ void CalcRhat() { #endif //KS: loop over chains and draws are independent so might as well collapse for sweet cache hits //Calculate the mean for each parameter within each considered chain + // MJR: Calculate using running totals to massively save on time and memory for (int m = 0; m < Nchains; ++m) { for (int j = 0; j < nDraw; ++j) From b25c5b007f5994fed9bc0335910661e56ddaa302 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 14:32:42 -0500 Subject: [PATCH 05/13] Updated data loading in original RHat macro to improve speed there in case people want to still use it --- Diagnostics/RHat_HighMem.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp index a0f9976d2..bc3910833 100644 --- a/Diagnostics/RHat_HighMem.cpp +++ b/Diagnostics/RHat_HighMem.cpp @@ -274,6 +274,14 @@ void PrepareChains() { } } + // MJR: array to hold branch values; SetBranchAddress in every step is very + // expensive, so doing it once only here saves time + double* branch_values = new double[nDraw](); + for (int j = 0; j < nDraw; ++j) + { + Chain->SetBranchAddress(BranchNames[j].Data(), &branch_values[j]); + } + //TN: move looping over toys here, so we don't need to loop over chains more than once if(BurnIn[m] >= nEntries[m]) { @@ -308,9 +316,8 @@ void PrepareChains() { // Set the branch addresses for params for (int j = 0; j < nDraw; ++j) { - Chain->SetBranchAddress(BranchNames[j].Data(), &Draws[m][i][j]); + Draws[m][i][j] = branch_values[j]; } - Chain->GetEntry(entry); }//end loop over toys From a99a293d5b83b5859bf9bda31e0808de32425454 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 15:10:48 -0500 Subject: [PATCH 06/13] Re-added thinning parameter to take place of the NToys command line argument --- Diagnostics/RHat.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 1e6d5398f..75c95166a 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -102,12 +102,14 @@ int main(int argc, char *argv[]) { if (argc < 2) { MACH3LOG_ERROR("Wrong arguments"); - MACH3LOG_ERROR("./RHat MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); + MACH3LOG_ERROR("./RHat NThin MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]"); throw MaCh3Exception(__FILE__ , __LINE__ ); } + NThin = atoi(argv[1]); + //KS Gelman suggests to diagnose on more than one chain - for (int i = 1; i < argc; i++) + for (int i = 2; i < argc; i++) { MCMCFile.push_back(std::string(argv[i])); MACH3LOG_INFO("Adding file: {}", MCMCFile.back()); @@ -166,7 +168,7 @@ void PrepareChains() { Chain->Add(MCMCFile[m].c_str()); nEntries[m] = int(Chain->GetEntries()); - Ntoys[m] = nEntries[m]; + Ntoys[m] = nEntries[m]/NThin; TotToys += Ntoys[m]; MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); @@ -283,7 +285,7 @@ void PrepareChains() { for (int i = 0; i < Ntoys[m]; i++) { // This is here as a placeholder in case we want to do some thinning later - int entry = i; + int entry = i*NThin; Chain->GetEntry(entry); From 88ffbc065e8d886881f71d9dff41484bbded692e Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 15:51:35 -0500 Subject: [PATCH 07/13] Added skeleton for FFT autocorrelation function support --- mcmc/MCMCProcessor.cpp | 7 ++++++- mcmc/MCMCProcessor.h | 6 ++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index 17a4ce8ba..d6181db9b 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -50,6 +50,7 @@ MCMCProcessor::MCMCProcessor(const std::string &InputFile) : ApplySmoothing = true; FancyPlotNames = true; doDiagMCMC = false; + useFFTAutoCorrelation = true; OutputSuffix = "_Process"; Post2DPlotThreshold = 1.e-5; @@ -3211,7 +3212,11 @@ void MCMCProcessor::DiagMCMC() { BatchedMeans(); // Draw the auto-correlations - AutoCorrelation(); + if (useFFTAutoCorrelation) { + AutoCorrelation_FFT(); + } else { + AutoCorrelation(); + } // Calculate Power Spectrum for each param PowerSpectrumAnalysis(); diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index ef91d7a00..5dec0fcdb 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -234,6 +234,8 @@ class MCMCProcessor { /// @brief Code will only plot 2D posteriors if Correlation are larger than defined threshold /// @param Threshold This threshold is compared with correlation value inline void SetPost2DPlotThreshold(const double Threshold){Post2DPlotThreshold = Threshold; }; + /// @brief Toggle using the FFT-based autocorrelation calculator + inline void SetUseFFTAutoCorrelation(const bool useFFT){useFFTAutoCorrelation = useFFT; }; /// @brief Setter related what parameters we want to exclude from analysis, for example if cross-section parameters look like xsec_, then passing "xsec_" will /// @param Batches Vector with parameters type names we want to exclude @@ -287,6 +289,8 @@ class MCMCProcessor { inline void ParamTraces(); /// @brief KS: Calculate autocorrelations supports both OpenMP and CUDA :) inline void AutoCorrelation(); + /// @brief MJR: Autocorrelation function using FFT algorithm for extra speed + inline void AutoCorrelation_FFT(); /// @brief KS: calc Effective Sample Size /// @param nLags Should be the same nLags as used in AutoCorrelation() /// @param LagL Value of LagL for each dial and each Lag @@ -405,6 +409,8 @@ class MCMCProcessor { bool ApplySmoothing; /// KS: Set Threshold when to plot 2D posterior as by default we get a LOT of plots double Post2DPlotThreshold; + /// MJR: Use FFT-based autocorrelation algorithm (save time & resources)? + bool useFFTAutoCorrelation; std::vector< int > NDSamplesBins; std::vector< std::string > NDSamplesNames; From 544cf9ea658c9f57e5370093274ec7f218e5287f Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 17:03:14 -0500 Subject: [PATCH 08/13] Implemented FFT-based autocorrelation calculator --- mcmc/MCMCProcessor.cpp | 105 +++++++++++++++++++++++++++++++++++++++++ mcmc/MCMCProcessor.h | 1 + 2 files changed, 106 insertions(+) diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index d6181db9b..800098622 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -3490,6 +3490,111 @@ void MCMCProcessor::ParamTraces() { OutputFile->cd(); } +// ********************************* +// MJR: Calculate autocorrelations using the FFT algorithm. +// Fast, even on CPU, and get all lags for free. +void MCMCProcessor::AutoCorrelation_FFT() { +// ********************************* + if (ParStep == nullptr) PrepareDiagMCMC(); + + TStopwatch clock; + clock.Start(); + const int nLags = AutoCorrLag; + MACH3LOG_INFO("Making auto-correlations for nLags = {}", nLags); + + // Prep outputs + OutputFile->cd(); + TDirectory* AutoCorrDir = OutputFile->mkdir("Auto_corr"); + std::vector LagKPlots(nDraw); + std::vector> LagL(nDraw); + + // Arrays needed to perform FFT using ROOT + double* ACFFT = new double[nEntries](); // Main autocorrelation array + double* ParVals = new double[nEntries](); // Param values for full chain + double* ParValsFFTR = new double[nEntries](); // FFT Real part + double* ParValsFFTI = new double[nEntries](); // FFT Imaginary part + double* ParValsFFTSquare = new double[nEntries](); // FFT Absolute square + double* ParValsComplex = new double[nEntries](); // Input Imaginary values (0) + + // Create forward and reverse FFT objects. I don't love using ROOT here, + // but it works so I can't complain + TVirtualFFT* fftf = TVirtualFFT::FFT(1, &nEntries, "C2CFORWARD"); + TVirtualFFT* fftb = TVirtualFFT::FFT(1, &nEntries, "C2CBACKWARD"); + + // Loop over all pars and calculate the full autocorrelation function using FFT + for (int j = 0; j < nDraw; ++j) { + // Initialize + LagL[j].resize(nLags); + for (int i = 0; i < nEntries; ++i) { + ParVals[i] = ParStep[j][i]-ParamSums[j]; // Subtract the mean to make it numerically tractable + ParValsComplex[i] = 0.; // Reset dummy array + } + + // Transform + fftf->SetPointsComplex(ParVals, ParValsComplex); + fftf->Transform(); + fftf->GetPointsComplex(ParValsFFTR, ParValsFFTI); + + // Square the results to get the power spectrum + for (int i = 0; i < nEntries; ++i) { + ParValsFFTSquare[i] = ParValsFFTR[i]*ParValsFFTR[i] + ParValsFFTI[i]*ParValsFFTI[i]; + } + + // Transforming back gives the autocovariance + fftb->SetPointsComplex(ParValsFFTSquare, ParValsComplex); + fftb->Transform(); + fftb->GetPointsComplex(ACFFT, ParValsComplex); + + // Divide by norm to get autocorrelation + double normAC = ACFFT[0]; + for (int i = 0; i < nEntries; ++ientry) { + ACFFT[i] /= normAC; + } + + // Get plotting info + TString Title = ""; + double Prior = 1.0; + double PriorError = 1.0; + GetNthParameter(j, Prior, PriorError, Title); + std::string HistName = Form("%s_%s_Lag", Title.Data(), BranchNames[j].Data()); + + // Initialize Lag plot + LagKPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nLags, 0.0, nLags); + LagKPlots[j]->GetXaxis()->SetTitle("Lag"); + LagKPlots[j]->GetYaxis()->SetTitle("Auto-correlation function"); + + // Fill plot + for (int k = 0; k < nLags; ++k) { + LagL[j][k] = ACFFT[k]; + LagKPlots[j]->SetBinContent(k, ACFFT[k]); + } + + // Write and clean up + AutoCorrDir>cd(); + LagKPlots[j]->Write(); + delete LagKPlots[j]; + } + + //KS: This is different diagnostic however it relies on calculated Lag, thus we call it before we delete LagKPlots + CalculateESS(nLags, LagL); + + // Clean up + delete[] ACFFT; + delete[] ParVals; + delete[] ParValsFFTR; + delete[] ParValsFFTI; + delete[] ParValsFFTSquare; + delete[] ParValsComplex; + + AutoCorrDir->Close(); + delete AutoCorrDir; + + OutputFile->cd(); + + clock.Stop(); + MACH3LOG_INFO("Making auto-correlations took {:.2f}s", clock.RealTime()); +} + // ********************************* //KS: Calculate autocorrelations supports both OpenMP and CUDA :) void MCMCProcessor::AutoCorrelation() { diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 5dec0fcdb..8fd43ed85 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -34,6 +34,7 @@ #include "TCandle.h" #include "TMath.h" #include "TMatrixDSymEigen.h" +#include "TVirtualFFT.h" // MaCh3 includes #include "mcmc/StatisticalUtils.h" From 5e081f03ed8e75d4130927593d3a45e2b248cbb6 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Fri, 22 Nov 2024 17:08:25 -0500 Subject: [PATCH 09/13] Fixed small compiler errors --- mcmc/MCMCProcessor.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index 800098622..a981100a0 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -3547,7 +3547,7 @@ void MCMCProcessor::AutoCorrelation_FFT() { // Divide by norm to get autocorrelation double normAC = ACFFT[0]; - for (int i = 0; i < nEntries; ++ientry) { + for (int i = 0; i < nEntries; ++i) { ACFFT[i] /= normAC; } @@ -3570,7 +3570,7 @@ void MCMCProcessor::AutoCorrelation_FFT() { } // Write and clean up - AutoCorrDir>cd(); + AutoCorrDir->cd(); LagKPlots[j]->Write(); delete LagKPlots[j]; } From 1b505c7e3076bd0f96af2dde3fb0b9b45502c2a8 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Mon, 25 Nov 2024 16:07:21 -0500 Subject: [PATCH 10/13] Forgot to clear branch_values array when done --- Diagnostics/RHat_HighMem.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp index bc3910833..47d2744f6 100644 --- a/Diagnostics/RHat_HighMem.cpp +++ b/Diagnostics/RHat_HighMem.cpp @@ -323,6 +323,7 @@ void PrepareChains() { //TN: There, we now don't need to keep the chain in memory anymore delete Chain; + delete[] branch_values; } //KS: Now prepare folded draws, quoting Gelman From 588964d595e91cafb2949b8c99bdc3ce443b91d1 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Mon, 25 Nov 2024 16:52:38 -0500 Subject: [PATCH 11/13] Changed RHat plots to not be fixed domain since stdev can vary wildly for parameters --- Diagnostics/RHat.cpp | 12 ++++++------ Diagnostics/RHat_HighMem.cpp | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 75c95166a..ee1a707a8 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -512,9 +512,9 @@ void SaveResults() { DiagFile->cd(); - TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", 200, 0, 2); - TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", 200, 0, 2); - TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", 200, 0, 2); + TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", nDraw, 0, nDraw); + TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", nDraw, 0, nDraw); + TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", nDraw, 0, nDraw); TH1D *RhatPlot = new TH1D("RhatPlot", "RhatPlot", 200, 0, 2); TH1D *EffectiveSampleSizePlot = new TH1D("EffectiveSampleSizePlot", "EffectiveSampleSizePlot", 400, 0, 10000); @@ -526,9 +526,9 @@ void SaveResults() { //KS: Fill only valid parameters if(ValidPar[j]) { - StandardDeviationGlobalPlot->Fill(StandardDeviationGlobal[j]); - BetweenChainVariancePlot->Fill(BetweenChainVariance[j]); - MarginalPosteriorVariancePlot->Fill(MarginalPosteriorVariance[j]); + StandardDeviationGlobalPlot->Fill(j,StandardDeviationGlobal[j]); + BetweenChainVariancePlot->Fill(j,BetweenChainVariance[j]); + MarginalPosteriorVariancePlot->Fill(j,MarginalPosteriorVariance[j]); RhatPlot->Fill(RHat[j]); EffectiveSampleSizePlot->Fill(EffectiveSampleSize[j]); if(RHat[j] > 1.1) Criterium++; diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp index 47d2744f6..fb67a56f3 100644 --- a/Diagnostics/RHat_HighMem.cpp +++ b/Diagnostics/RHat_HighMem.cpp @@ -614,15 +614,15 @@ void SaveResults() { DiagFile->cd(); - TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", 200, 0, 2); - TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", 200, 0, 2); - TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", 200, 0, 2); + TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", nDraw, 0, nDraw); + TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", nDraw, 0, nDraw); + TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", nDraw, 0, nDraw); TH1D *RhatPlot = new TH1D("RhatPlot", "RhatPlot", 200, 0, 2); TH1D *EffectiveSampleSizePlot = new TH1D("EffectiveSampleSizePlot", "EffectiveSampleSizePlot", 400, 0, 10000); - TH1D *StandardDeviationGlobalFoldedPlot = new TH1D("StandardDeviationGlobalFoldedPlot", "StandardDeviationGlobalFoldedPlot", 200, 0, 2); - TH1D *BetweenChainVarianceFoldedPlot = new TH1D("BetweenChainVarianceFoldedPlot", "BetweenChainVarianceFoldedPlot", 200, 0, 2); - TH1D *MarginalPosteriorVarianceFoldedPlot = new TH1D("MarginalPosteriorVarianceFoldedPlot", "MarginalPosteriorVarianceFoldedPlot", 200, 0, 2); + TH1D *StandardDeviationGlobalFoldedPlot = new TH1D("StandardDeviationGlobalFoldedPlot", "StandardDeviationGlobalFoldedPlot", nDraw, 0, nDraw); + TH1D *BetweenChainVarianceFoldedPlot = new TH1D("BetweenChainVarianceFoldedPlot", "BetweenChainVarianceFoldedPlot", nDraw, 0, nDraw); + TH1D *MarginalPosteriorVarianceFoldedPlot = new TH1D("MarginalPosteriorVarianceFoldedPlot", "MarginalPosteriorVarianceFoldedPlot", nDraw, 0, nDraw); TH1D *RhatFoldedPlot = new TH1D("RhatFoldedPlot", "RhatFoldedPlot", 200, 0, 2); TH1D *EffectiveSampleSizeFoldedPlot = new TH1D("EffectiveSampleSizeFoldedPlot", "EffectiveSampleSizeFoldedPlot", 400, 0, 10000); From f8c1fbb4018150057721c1ebb4be66c8abc4ed15 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Mon, 25 Nov 2024 17:41:01 -0500 Subject: [PATCH 12/13] Fixed a bug where I was accidentally overwriting the number of generated toys at each step, whoops --- Diagnostics/RHat.cpp | 36 +++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index ee1a707a8..2a45f9997 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -40,7 +40,8 @@ /// // ******************* -int* Ntoys; +int* Ntoys_requested; +int* Ntoys_filled; int TotToys; int NThin; int Nchains; @@ -145,7 +146,8 @@ void PrepareChains() { TStopwatch clock; clock.Start(); - Ntoys = new int[Nchains](); + Ntoys_requested = new int[Nchains](); + Ntoys_filled = new int[Nchains](); TotToys = 0; std::vector BurnIn(Nchains); std::vector nEntries(Nchains); @@ -168,11 +170,11 @@ void PrepareChains() { Chain->Add(MCMCFile[m].c_str()); nEntries[m] = int(Chain->GetEntries()); - Ntoys[m] = nEntries[m]/NThin; - TotToys += Ntoys[m]; + Ntoys_requested[m] = nEntries[m]/NThin; + Ntoys_filled[m] = 0; MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str()); - MACH3LOG_INFO("Generating {}", Ntoys[m]); + MACH3LOG_INFO("Generating {} Toys", Ntoys_requested[m]); // Set the step cut to be 20% BurnIn[m] = nEntries[m]/5; @@ -282,7 +284,7 @@ void PrepareChains() { } MACH3LOG_INFO("Loading chain {} / {}...", m, Nchains); - for (int i = 0; i < Ntoys[m]; i++) + for (int i = 0; i < Ntoys_requested[m]; i++) { // This is here as a placeholder in case we want to do some thinning later int entry = i*NThin; @@ -293,13 +295,12 @@ void PrepareChains() { // Note, entry is not necessarily the same as the step due to merged ROOT files, so can't choose an entry in the range BurnIn - nEntries :( if (step[m] < BurnIn[m]) { - Ntoys[m]--; continue; } // Output some info for the user - if (Ntoys[m] > 10 && i % (Ntoys[m]/10) == 0) { - MaCh3Utils::PrintProgressBar(i+m*Ntoys[m], static_cast(Ntoys[m])*Nchains); + if (Ntoys_requested[m] > 10 && i % (Ntoys_requested[m]/10) == 0) { + MaCh3Utils::PrintProgressBar(i+m*Ntoys_requested[m], static_cast(Ntoys_requested[m])*Nchains); MACH3LOG_DEBUG("Getting random entry {}", entry); } @@ -314,6 +315,10 @@ void PrepareChains() { S2_chain[m][j] += branch_values[j]*branch_values[j]; } + // Increment counters + Ntoys_filled[m]++; + TotToys++; + }//end loop over toys //TN: There, we now don't need to keep the chain in memory anymore @@ -399,8 +404,8 @@ void CalcRhat() { { for (int j = 0; j < nDraw; ++j) { - Mean[m][j] += S1_chain[m][j] / static_cast(Ntoys[m]); - StandardDeviation[m][j] = S2_chain[m][j]/static_cast(Ntoys[m]) - Mean[m][j]*Mean[m][j]; + Mean[m][j] = S1_chain[m][j] / static_cast(Ntoys_filled[m]); + StandardDeviation[m][j] = S2_chain[m][j]/static_cast(Ntoys_filled[m]) - Mean[m][j]*Mean[m][j]; } } @@ -412,10 +417,6 @@ void CalcRhat() { { for (int m = 0; m < Nchains; ++m) { - if (m == 0) - { - StandardDeviationGlobal[j] = 0.0; - } StandardDeviationGlobal[j] += StandardDeviation[m][j]; } MeanGlobal[j] = S1_global[j] / static_cast(TotToys); @@ -436,7 +437,7 @@ void CalcRhat() { { for (int m = 0; m < Nchains; ++m) { - BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]) * Ntoys[m]; + BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]) * Ntoys_filled[m]; } BetweenChainVariance[j] = BetweenChainVariance[j]/(Nchains-1); } @@ -671,7 +672,8 @@ void DestroyArrays() { delete[] S1_global; delete[] S2_global; - delete[] Ntoys; + delete[] Ntoys_requested; + delete[] Ntoys_filled; } // ******************* From ad609396686696186906e7085a8cf87e896b47d2 Mon Sep 17 00:00:00 2001 From: mreh528 Date: Mon, 25 Nov 2024 17:43:49 -0500 Subject: [PATCH 13/13] Forgot to fix plotting for HighMem RHat --- Diagnostics/RHat_HighMem.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp index fb67a56f3..575ae900d 100644 --- a/Diagnostics/RHat_HighMem.cpp +++ b/Diagnostics/RHat_HighMem.cpp @@ -636,17 +636,17 @@ void SaveResults() { //KS: Fill only valid parameters if(ValidPar[j]) { - StandardDeviationGlobalPlot->Fill(StandardDeviationGlobal[j]); - BetweenChainVariancePlot->Fill(BetweenChainVariance[j]); - MarginalPosteriorVariancePlot->Fill(MarginalPosteriorVariance[j]); + StandardDeviationGlobalPlot->Fill(j,StandardDeviationGlobal[j]); + BetweenChainVariancePlot->Fill(j,BetweenChainVariance[j]); + MarginalPosteriorVariancePlot->Fill(j,MarginalPosteriorVariance[j]); RhatPlot->Fill(RHat[j]); EffectiveSampleSizePlot->Fill(EffectiveSampleSize[j]); if(RHat[j] > 1.1) Criterium++; - StandardDeviationGlobalFoldedPlot->Fill(StandardDeviationGlobalFolded[j]); - BetweenChainVarianceFoldedPlot->Fill(BetweenChainVarianceFolded[j]); - MarginalPosteriorVarianceFoldedPlot->Fill(MarginalPosteriorVarianceFolded[j]); + StandardDeviationGlobalFoldedPlot->Fill(j,StandardDeviationGlobalFolded[j]); + BetweenChainVarianceFoldedPlot->Fill(j,BetweenChainVarianceFolded[j]); + MarginalPosteriorVarianceFoldedPlot->Fill(j,MarginalPosteriorVarianceFolded[j]); RhatFoldedPlot->Fill(RHatFolded[j]); EffectiveSampleSizeFoldedPlot->Fill(EffectiveSampleSizeFolded[j]); if(RHatFolded[j] > 1.1) CiteriumFolded++;