Merge pull request #102 from mach3-software/feature_gif

Paramter Evolution .gif

Diagnostics/Diag_Config.yaml

@@ -58,6 +58,7 @@ ProcessMCMC:
   CalcBayesFactor: true
   CalcSavageDickey: true
   CalcBipolarPlot: false
+  CalcParameterEvolution: false
 
   #Functionality 2D
   PlotCorr: false
@@ -119,6 +120,10 @@ ProcessMCMC:
     - ["CCQE_Eb", ["EB_dial_C_nu", "EB_dial_O_nu", "EB_dial_C_nubar", "EB_dial_O_nubar", "Alpha_q3"]]
     - ["Oscillation", ["delta_cp", "delm2_23", "sin2th_13", "sin2th_23"]]
 
+  #First is paramter name and second number of intervals
+  ParameterEvolution:
+    - ["EB_dial_C_nu", 20]
+
   #First param name, then new prior and new prior error
   # apply reactor constraints [vector of param names] [vector of new central value] [vector of new prior errors]
   PriorReweighting: [ ["sin2th_13"], [0.0220], [0.0007] ]

Diagnostics/ProcessMCMC.cpp

@@ -9,6 +9,7 @@ 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);
@@ -105,6 +106,7 @@ void ProcessMCMC(const std::string& inputFile)
   if(GetFromManager<bool>(Settings["CalcBayesFactor"], true))  CalcBayesFactor(Processor);
   if(GetFromManager<bool>(Settings["CalcSavageDickey"], true)) CalcSavageDickey(Processor);
   if(GetFromManager<bool>(Settings["CalcBipolarPlot"], false)) CalcBipolarPlot(Processor);
+  if(GetFromManager<bool>(Settings["CalcParameterEvolution"], false)) CalcParameterEvolution(Processor);
 
   if(PlotCorr)
   {
@@ -361,6 +363,23 @@ void CalcSavageDickey(MCMCProcessor* Processor)
   return;
 }
 
+void CalcParameterEvolution(MCMCProcessor* Processor)
+{
+  YAML::Node card_yaml = YAML::LoadFile(config.c_str());
+  YAML::Node Settings = card_yaml["ProcessMCMC"];
+
+  std::vector<std::string> ParNames;
+  std::vector<int> Intervals;
+
+  for (const auto& d : Settings["ParameterEvolution"])
+  {
+    ParNames.push_back(d[0].as<std::string>());
+    Intervals.push_back(d[1].as<int>());
+  }
+  Processor->ParameterEvolution(ParNames, Intervals);
+  return;
+}
+
 void CalcBipolarPlot(MCMCProcessor* Processor)
 {
   YAML::Node card_yaml = YAML::LoadFile(config.c_str());
@@ -376,6 +395,7 @@ void CalcBipolarPlot(MCMCProcessor* Processor)
   return;
 }
 
+
 void GetTrianglePlot(MCMCProcessor* Processor) {
   YAML::Node card_yaml = YAML::LoadFile(config.c_str());
   YAML::Node Settings = card_yaml["ProcessMCMC"];

mcmc/MCMCProcessor.cpp

@@ -3195,6 +3195,76 @@ void MCMCProcessor::ReweightPrior(const std::vector<std::string>& Names,
   OutputFile->cd();
 }
 
+// **************************
+// Diagnose the MCMC
+void MCMCProcessor::ParameterEvolution(const std::vector<std::string>& Names,
+                                       const std::vector<int>& NIntervals) {
+// **************************
+  MACH3LOG_INFO("Parameter Evolution gif");
+
+  //KS: First we need to find parameter number based on name
+  for(unsigned int k = 0; k < Names.size(); ++k)
+  {
+    //KS: First we need to find parameter number based on name
+    int ParamNo = GetParamIndexFromName(Names[k]);
+    if(ParamNo == _UNDEF_)
+    {
+      MACH3LOG_WARN("Couldn't find param {}. Can't reweight Prior", Names[k]);
+      continue;
+    }
+
+    const int IntervalsSize = nSteps/NIntervals[k];
+
+    // ROOT won't overwrite gifs so we need to delete the file if it's there already
+    std::remove(std::string(Names[k]+".gif").c_str());
+
+    // This holds the posterior density
+    const double maxi = Chain->GetMaximum(BranchNames[ParamNo]);
+    const double mini = Chain->GetMinimum(BranchNames[ParamNo]);
+
+    int Counter = 0;
+    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);
+      EvePlot->SetMinimum(0);
+      EvePlot->GetYaxis()->SetTitle("PDF");
+      EvePlot->GetYaxis()->SetNoExponent(false);
+
+      //KS: Apply additional Cuts, like mass ordering
+      std::string CutPosterior1D = "step > " + std::to_string(i*IntervalsSize+IntervalsSize);
+
+      std::string TextTitle = "Steps = 0 - "+std::to_string(Counter*IntervalsSize+IntervalsSize);
+      // Project BranchNames[ParamNo] onto hpost, applying stepcut
+      Chain->Project(BranchNames[ParamNo], BranchNames[ParamNo], CutPosterior1D.c_str());
+
+      EvePlot->SetLineWidth(2);
+      EvePlot->SetLineColor(kBlue-1);
+      EvePlot->SetTitle(Names[k].c_str());
+      EvePlot->GetXaxis()->SetTitle(EvePlot->GetTitle());
+      EvePlot->GetYaxis()->SetLabelOffset(1000);
+      if(ApplySmoothing) EvePlot->Smooth();
+
+      EvePlot->Scale(1. / EvePlot->Integral());
+
+      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");
+
+      if(i == 0) Posterior->Print((std::string(Names[k] + ".gif++20").c_str())); // produces infinite loop animated GIF
+      else Posterior->Print((std::string(Names[k]+".gif+20").c_str())); // add picture to .gif
+
+      delete EvePlot;
+      delete text;
+      Counter++;
+    }
+  }
+}
+
 // **************************
 // Diagnose the MCMC
 void MCMCProcessor::DiagMCMC() {

mcmc/MCMCProcessor.h

@@ -6,6 +6,7 @@
 
 // C++ includes
 #include <complex>
+#include <cstdio>
 
 // ROOT includes
 #include "TObjArray.h"
@@ -152,6 +153,12 @@ class MCMCProcessor {
                        const std::vector<double>& NewCentral,
                        const std::vector<double>& NewError);
 
+    /// @brief Make .gif of parameter evolution
+    /// @param ParName Parameter names for which we do .gif
+    /// @param NIntervals Number of intervals for a gif
+    void ParameterEvolution(const std::vector<std::string>& Names,
+                            const std::vector<int>& NIntervals);
+
     /// @brief KS: Perform MCMC diagnostic including Autocorrelation, Trace etc.
     void DiagMCMC();
 

mcmc/gpuMCMCProcessorUtils.cuh

@@ -9,10 +9,6 @@
 
 #include "manager/gpuUtils.cuh"
 
-// *******************************************
-//              INITIALISE GPU
-// *******************************************
-
 /// @brief KS: Initialiser, here we allocate memory for variables and copy constants
 /// @param ParStep_gpu Parameter value at each step
 /// @param NumeratorSum_gpu Sum used for nominator of autocorrelation calculations