Implementation of Prior/Posterior Predictive and Bayesian p-Value calculations following the approach described in [10]. More...

#include <Fitters/PredictiveThrower.h>

Inheritance diagram for PredictiveThrower:

[legend]

Collaboration diagram for PredictiveThrower:

[legend]

Public Member Functions
	PredictiveThrower (manager *const fitMan)
	Constructor. More...

virtual	~PredictiveThrower ()
	Destructor. More...

void	ProduceToys ()
	Produce toys by throwing from MCMC. More...

void	RunPredictiveAnalysis ()
	Main routine responsible for producing posterior predictive distributions and $p$-value. More...

void	RunMCMC () override
	This is not used in this class. More...

Public Member Functions inherited from FitterBase
	FitterBase (manager *const fitMan)
	Constructor. More...

virtual	~FitterBase ()
	Destructor for the FitterBase class. More...

void	AddSampleHandler (SampleHandlerBase *sample)
	This function adds a sample PDF object to the analysis framework. The sample PDF object will be utilized in fitting procedures or likelihood scans. More...

void	AddSystObj (ParameterHandlerBase *cov)
	This function adds a Covariance object to the analysis framework. The Covariance object will be utilized in fitting procedures or likelihood scans. More...

void	DragRace (const int NLaps=100)
	Calculates the required time for each sample or covariance object in a drag race simulation. Inspired by Dan's feature. More...

void	RunLLHScan ()
	Perform a 1D likelihood scan. More...

void	GetStepScaleBasedOnLLHScan ()
	LLH scan is good first estimate of step scale. More...

void	Run2DLLHScan ()
	Perform a 2D likelihood scan. More...

void	RunSigmaVar ()
	Perform a 2D and 1D sigma var for all samples. More...

void	RunSigmaVarFD ()
	Perform a 1D sigma var for all samples. More...

virtual void	StartFromPreviousFit (const std::string &FitName)
	Allow to start from previous fit/chain. More...

std::string	GetName () const
	Get name of class. More...

Private Member Functions
void	SetParamters ()
	This set some params to prior value this way you can evaluate errors from subset of errors. More...

void	SetupToyGeneration ()
	Setup useful variables etc before stating toy generation. More...

bool	LoadToys ()
	Load existing toys. More...

void	SetupSampleInformation ()
	Setup sample information. More...

std::unique_ptr< TH1D >	MakePredictive (const std::vector< std::unique_ptr< TH1D >> &Toys, const std::string &Sample_Name, const std::string &suffix, const bool DebugHistograms)
	Produce posterior predictive distribution. More...

std::vector< std::unique_ptr< TH2D > >	ProduceSpectra (const std::vector< std::vector< std::unique_ptr< TH1D >>> &Toys, const std::string suffix)
	Produce Violin style spectra. More...

void	PosteriorPredictivepValue (const std::vector< std::unique_ptr< TH1D >> &PostPred_mc, const std::vector< TDirectory * > &SampleDir)
	Calculate Posterior Predictive $p$-value. More...

double	GetLLH (const std::unique_ptr< TH1D > &DatHist, const std::unique_ptr< TH1D > &MCHist, const std::unique_ptr< TH1D > &W2Hist, SampleHandlerBase *SampleHandler)
	Helper functions to calculate likelihoods using TH1D. More...

void	MakeChi2Plots (const std::vector< std::vector< double >> &Chi2_x, const std::string &Chi2_x_title, const std::vector< std::vector< double >> &Chi2_y, const std::string &Chi2_y_title, const std::vector< TDirectory * > &SampleDir, const std::string Title)
	Produce Chi2 plot for a single sample based on which $p$-value is calculated. More...

Private Attributes
bool	FullLLH
	KS: Use Full LLH or only sample contribution based on discussion with Asher we almost always only want the sample likelihood. More...

int	NModelParams
	KS: Count total number of model parameters which can be used for stuff like BIC. More...

bool	Is_PriorPredictive
	Whether it is Prior or Posterior predictive. More...

int	TotalNumberOfSamples
	Number of toys we are generating analysing. More...

std::vector< std::string >	SampleNames
	Name of a single sample. More...

std::vector< int >	SampleObjectMap
	Maps if sample with given SampleHandler, useful if we have more than one sample in single object. More...

int	Ntoys
	Number of toys we are generating analysing. More...

std::vector< std::string >	ParameterGroupsNotVaried
	KS: Names of parameter groups that will not be varied. More...

std::unordered_set< int >	ParameterOnlyToVary
	KS: Index of parameters groups that will be varied. More...

ParameterHandlerGeneric *	ModelSystematic
	Pointer to El Generico. More...

std::vector< std::unique_ptr< TH1D > >	Data_Hist
	Vector of Data histograms. More...

std::vector< std::unique_ptr< TH1D > >	MC_Nom_Hist
	Vector of MC histograms. More...

std::vector< std::unique_ptr< TH1D > >	W2_Nom_Hist
	Vector of W2 histograms. More...

std::vector< std::vector< std::unique_ptr< TH1D > > >	MC_Hist_Toy

std::vector< std::vector< std::unique_ptr< TH1D > > >	W2_Hist_Toy

std::vector< double >	ReweightWeight
	Reweighting factors applied for each toy, by default 1. More...

std::vector< double >	PenaltyTerm
	Penalty term values for each toy by default 0. More...

Additional Inherited Members
Protected Member Functions inherited from FitterBase
void	ProcessMCMC ()
	Process MCMC output. More...

void	PrepareOutput ()
	Prepare the output file. More...

void	SaveOutput ()
	Save output and close files. More...

void	SanitiseInputs ()
	Remove obsolete memory and make other checks before fit starts. More...

void	SaveSettings ()
	Save the settings that the MCMC was run with. More...

bool	GetScaneRange (std::map< std::string, std::vector< double >> &scanRanges)
	YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow. More...

bool	CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
	KS: Check whether we want to skip parameter using skip vector. More...

void	CustomRange (const std::string &ParName, const double sigma, double &ParamShiftValue)
	For comparison with P-Theta we usually have to apply different parameter values then usual 1, 3 sigma. More...

Protected Attributes inherited from FitterBase
manager *	fitMan
	The manager. More...

unsigned int	step
	current state More...

double	logLCurr
	current likelihood More...

double	logLProp
	proposed likelihood More...

double	accProb
	current acceptance prob More...

int	accCount
	counts accepted steps More...

unsigned int	stepStart
	step start, by default 0 if we start from previous chain then it will be different More...

std::vector< double >	sample_llh
	store the llh breakdowns More...

std::vector< double >	syst_llh
	systematic llh breakdowns More...

std::vector< SampleHandlerBase * >	samples
	Sample holder. More...

unsigned int	TotalNSamples
	Total number of samples used. More...

std::vector< ParameterHandlerBase * >	systematics
	Systematic holder. More...

std::unique_ptr< TStopwatch >	clock
	tells global time how long fit took More...

std::unique_ptr< TStopwatch >	stepClock
	tells how long single step/fit iteration took More...

double	stepTime
	Time of single step. More...

std::unique_ptr< TRandom3 >	random
	Random number. More...

TFile *	outputFile
	Output. More...

TDirectory *	CovFolder
	Output cov folder. More...

TDirectory *	SampleFolder
	Output sample folder. More...

TTree *	outTree
	Output tree with posteriors. More...

int	auto_save
	auto save every N steps More...

bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO. More...

bool	FileSaved
	Checks if file saved not repeat some operations. More...

bool	SettingsSaved
	Checks if setting saved not repeat some operations. More...

bool	OutputPrepared
	Checks if output prepared not repeat some operations. More...

std::string	AlgorithmName
	Name of fitting algorithm that is being used. More...

Detailed Description

Implementation of Prior/Posterior Predictive and Bayesian p-Value calculations following the approach described in [10].

For more information, visit the Wiki.

Author: Asher Kaboth; Dan Barrow; Ed Atkin; Yashwanth S Prabhu; Kamil Skwarczynski; Patrick Dunne; Clarence Wret

Todo:

add BIC, DIC, WAIC

add ability yo make projection for Get1DDiscVar

add ability for TH2D

speed improvements

add Rate $p$-value

unify code with SampleSummary

Definition at line 23 of file PredictiveThrower.h.

Constructor & Destructor Documentation

◆ PredictiveThrower()

PredictiveThrower::PredictiveThrower ( manager *const fitMan )

Constructor.

Parameters

fitMan A pointer to a manager object, which will handle all settings.

Definition at line 6 of file PredictiveThrower.cpp.

                                                  : FitterBase(man) {
 // *************************
  AlgorithmName = "PredictiveThrower";
   if(!CheckNodeExists(fitMan->raw(), "Predictive")) {
    MACH3LOG_ERROR("Predictive is missing in your main yaml config");
    throw MaCh3Exception(__FILE__ , __LINE__ );
   }
  
   ModelSystematic = nullptr;
   // Use the full likelihood for the Prior/Posterior predictive pvalue
   FullLLH = GetFromManager<bool>(fitMan->raw()["Predictive"]["FullLLH"], false, __FILE__, __LINE__ );
   NModelParams = 0;
  
   Is_PriorPredictive = Get<bool>(fitMan->raw()["Predictive"]["PriorPredictive"], __FILE__, __LINE__);
   Ntoys = Get<int>(fitMan->raw()["Predictive"]["Ntoy"], __FILE__, __LINE__);
  
   ReweightWeight.resize(Ntoys);
   PenaltyTerm.resize(Ntoys);
 }

◆ ~PredictiveThrower()

PredictiveThrower::~PredictiveThrower ( )

virtual

Destructor.

Definition at line 28 of file PredictiveThrower.cpp.

                                       {
 // *************************
  
 }

Member Function Documentation

◆ GetLLH()

double PredictiveThrower::GetLLH	(	const std::unique_ptr< TH1D > &	DatHist,
		const std::unique_ptr< TH1D > &	MCHist,
		const std::unique_ptr< TH1D > &	W2Hist,
		SampleHandlerBase *	SampleHandler
	)

private

Helper functions to calculate likelihoods using TH1D.

Parameters

Data	histogram with data distribution for a single sample
MC	histogram with MC distribution for a single sample
W2	histogram with W2 distribution for a single sample

Definition at line 583 of file PredictiveThrower.cpp.

                                                                    {
 // *************************
   double llh = 0.0;
   for (int i = 1; i <= DatHist->GetXaxis()->GetNbins(); ++i)
   {
     const double data = DatHist->GetBinContent(i);
     const double mc = MCHist->GetBinContent(i);
     const double w2 = W2Hist->GetBinContent(i);
     llh += SampleHandler->GetTestStatLLH(data, mc, w2);
   }
   //KS: do times 2 because banff reports chi2
   return 2*llh;
 }

◆ LoadToys()

bool PredictiveThrower::LoadToys ( )

private

Load existing toys.

Definition at line 191 of file PredictiveThrower.cpp.

                                  {
 // *************************
   auto PosteriorFileName = Get<std::string>(fitMan->raw()["Predictive"]["PosteriorFile"], __FILE__, __LINE__);
   // Open the ROOT file
   TFile* file = TFile::Open(PosteriorFileName.c_str(), "READ");
   TDirectory* ToyDir = nullptr;
   if (!file || file->IsZombie()) {
     return false;
   } else {
     // Check for the "toys" directory
     if ((ToyDir = file->GetDirectory("Toys"))) {
       MACH3LOG_INFO("Found toys in Posterior file will attempt toy reading");
     } else {
       file->Close();
       delete file;
       return false;
     }
   }
  
   // Finally get the TTree branch with the penalty vectors for each of the toy throws
   TTree* PenaltyTree = static_cast<TTree*>(file->Get("ToySummary"));
   if (!PenaltyTree) {
     MACH3LOG_WARN("ToySummary TTree not found in file.");
     file->Close();
     delete file;
     return false;
   }
  
   Ntoys = static_cast<int>(PenaltyTree->GetEntries());
   int ConfigNtoys = Get<int>(fitMan->raw()["Predictive"]["Ntoy"], __FILE__, __LINE__);;
   if (Ntoys != ConfigNtoys) {
     MACH3LOG_WARN("Found different number of toys in saved file than asked to run!");
     MACH3LOG_INFO("Will read _ALL_ toys in the file");
     MACH3LOG_INFO("Ntoys in file: {}", Ntoys);
     MACH3LOG_INFO("Ntoys specified: {}", ConfigNtoys);
   }
  
   PenaltyTerm.resize(Ntoys);
   ReweightWeight.resize(Ntoys);
  
   double Penalty = 0, Weight = 1;
   PenaltyTree->SetBranchAddress("Penalty", &Penalty);
   PenaltyTree->SetBranchAddress("Weight", &Weight);
   PenaltyTree->SetBranchAddress("NModelParams", &NModelParams);
  
   for (int i = 0; i < Ntoys; ++i) {
     PenaltyTree->GetEntry(i);
     if (FullLLH) {
       PenaltyTerm[i] = Penalty;
     } else {
       PenaltyTerm[i] = 0.0;
     }
  
     ReweightWeight[i] = Weight;
   }
   // Resize all vectors and get sample names
   SetupSampleInformation();
  
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     TH1D* DataHist1D = static_cast<TH1D*>(ToyDir->Get((SampleNames[sample] + "_data").c_str()));
     Data_Hist[sample] = M3::Clone(DataHist1D);
  
     TH1D* MCHist1D = static_cast<TH1D*>(ToyDir->Get((SampleNames[sample] + "_mc").c_str()));
     MC_Nom_Hist[sample] = M3::Clone(MCHist1D);
  
     TH1D* W2Hist1D = static_cast<TH1D*>(ToyDir->Get((SampleNames[sample] + "_w2").c_str()));
     W2_Nom_Hist[sample] = M3::Clone(W2Hist1D);
   }
  
  
   for (int iToy = 0; iToy < Ntoys; ++iToy)
   {
     if (iToy % 100 == 0) MACH3LOG_INFO("   Loaded toy {}", iToy);
  
     for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
       TH1D* MCHist1D = static_cast<TH1D*>(ToyDir->Get((SampleNames[sample] + "_mc_" + std::to_string(iToy)).c_str()));
       TH1D* W2Hist1D = static_cast<TH1D*>(ToyDir->Get((SampleNames[sample] + "_w2_" + std::to_string(iToy)).c_str()));
  
       MC_Hist_Toy[sample][iToy] = M3::Clone(MCHist1D);
       W2_Hist_Toy[sample][iToy] = M3::Clone(W2Hist1D);
     }
   }
  
   file->Close();
   delete file;
   return true;
 }

◆ MakeChi2Plots()

void PredictiveThrower::MakeChi2Plots	(	const std::vector< std::vector< double >> &	Chi2_x,
		const std::string &	Chi2_x_title,
		const std::vector< std::vector< double >> &	Chi2_y,
		const std::string &	Chi2_y_title,
		const std::vector< TDirectory * > &	SampleDir,
		const std::string	Title
	)

private

Produce Chi2 plot for a single sample based on which $p$-value is calculated.

Definition at line 644 of file PredictiveThrower.cpp.

                                           {
 // *************************
   for (int iSample = 0; iSample < TotalNumberOfSamples+1; ++iSample) {
     SampleDir[iSample]->cd();
  
     // Transpose to extract chi2 values for a given sample across all toys
     std::vector<double> chi2_y_sample(Ntoys);
     std::vector<double> chi2_x_per_sample(Ntoys);
  
     for (int iToy = 0; iToy < Ntoys; ++iToy) {
       chi2_y_sample[iToy] = Chi2_y[iToy][iSample];
       chi2_x_per_sample[iToy]  = Chi2_x[iToy][iSample];
     }
  
     const double min_val = std::min(*std::min_element(chi2_y_sample.begin(), chi2_y_sample.end()),
                               *std::min_element(chi2_x_per_sample.begin(), chi2_x_per_sample.end()));
     const double max_val = std::max(*std::max_element(chi2_y_sample.begin(), chi2_y_sample.end()),
                               *std::max_element(chi2_x_per_sample.begin(), chi2_x_per_sample.end()));
  
     auto chi2_hist = std::make_unique<TH2D>((SampleNames[iSample] + Title).c_str(),
                                             (SampleNames[iSample] + Title).c_str(),
                                             100, min_val, max_val, 100, min_val, max_val);
     chi2_hist->SetDirectory(nullptr);
     chi2_hist->GetXaxis()->SetTitle(Chi2_x_title.c_str());
     chi2_hist->GetYaxis()->SetTitle(Chi2_y_title.c_str());
  
     for (int iToy = 0; iToy < Ntoys; ++iToy) {
       chi2_hist->Fill(chi2_x_per_sample[iToy], chi2_y_sample[iToy]);
     }
  
     Get2DBayesianpValue(chi2_hist.get());
     chi2_hist->Write();
   }
 }

◆ MakePredictive()

std::unique_ptr< TH1D > PredictiveThrower::MakePredictive	(	const std::vector< std::unique_ptr< TH1D >> &	Toys,
		const std::string &	Sample_Name,
		const std::string &	suffix,
		const bool	DebugHistograms
	)

private

Produce posterior predictive distribution.

Definition at line 494 of file PredictiveThrower.cpp.

                                                                                     {
 // *************************
   constexpr int nXBins = 100;
   int nbinsx = Toys[0]->GetNbinsX();
  
   auto PredictiveHist = std::make_unique<TH1D>((Sample_Name + "_" + suffix + "_PostPred").c_str(),
                                                (Sample_Name + "_" + suffix + "_PostPred").c_str(),
                                                 nbinsx, Toys[0]->GetXaxis()->GetXbins()->GetArray());
   PredictiveHist->SetDirectory(nullptr);
  
   for (int i = 1; i <= nbinsx; ++i) {
     double x_low  = Toys[0]->GetXaxis()->GetBinLowEdge(i);
     double x_high = Toys[0]->GetXaxis()->GetBinUpEdge(i);
     TString projName = TString::Format("%s %s X: [%.3f, %.3f]", Sample_Name.c_str(), suffix.c_str(), x_low, x_high);
     auto PosteriorHist = std::make_unique<TH1D>(projName, projName, nXBins, 1, -1);
     PosteriorHist->SetDirectory(nullptr);
  
     for (size_t iToy = 0; iToy < Toys.size(); ++iToy) {
       const double Content = Toys[iToy]->GetBinContent(i);
       PosteriorHist->Fill(Content, ReweightWeight[iToy]);
     }
  
     if(DebugHistograms) PosteriorHist->Write();
  
     const double nMean = PosteriorHist->GetMean();
     const double nMeanError = PosteriorHist->GetRMS();
  
     PredictiveHist->SetBinContent(i, nMean);
     PredictiveHist->SetBinError(i, nMeanError);
   }
  
   PredictiveHist->Write();
   return PredictiveHist;
 }

◆ PosteriorPredictivepValue()

void PredictiveThrower::PosteriorPredictivepValue	(	const std::vector< std::unique_ptr< TH1D >> &	PostPred_mc,
		const std::vector< TDirectory * > &	SampleDir
	)

private

Calculate Posterior Predictive $p$-value.

TODO This can be multithreaded

Definition at line 601 of file PredictiveThrower.cpp.

                                                                                            {
 // *************************
   //(void) PostPred_w2;
   // [Toys][Sample]
   std::vector<std::vector<double>> chi2_dat_vec(Ntoys);
   std::vector<std::vector<double>> chi2_mc_vec(Ntoys);
   std::vector<std::vector<double>> chi2_pred_vec(Ntoys);
  
   for(int iToy = 0; iToy < Ntoys; iToy++) {
     chi2_dat_vec[iToy].resize(TotalNumberOfSamples+1, 0);
     chi2_mc_vec[iToy].resize(TotalNumberOfSamples+1, 0);
     chi2_pred_vec[iToy].resize(TotalNumberOfSamples+1, 0);
  
     chi2_dat_vec[iToy].back() = PenaltyTerm[iToy];
     chi2_mc_vec[iToy].back() = PenaltyTerm[iToy];
     chi2_pred_vec[iToy].back() = PenaltyTerm[iToy];
  
     for (int iSample = 0; iSample < TotalNumberOfSamples; ++iSample) {
       auto DrawFluctHist = M3::Clone(MC_Hist_Toy[iSample][iToy].get());
       MakeFluctuatedHistogramAlternative(DrawFluctHist.get(), MC_Hist_Toy[iSample][iToy].get(), random.get());
  
       auto PredFluctHist = M3::Clone(PostPred_mc[iSample].get());
       MakeFluctuatedHistogramAlternative(PredFluctHist.get(), PostPred_mc[iSample].get(), random.get());
  
       // Okay now we can do our chi2 calculation for our sample
       chi2_dat_vec[iToy][iSample]  = GetLLH(Data_Hist[iSample], MC_Hist_Toy[iSample][iToy], W2_Hist_Toy[iSample][iToy], samples[SampleObjectMap[iSample]]);
       chi2_mc_vec[iToy][iSample]   = GetLLH(DrawFluctHist, MC_Hist_Toy[iSample][iToy], W2_Hist_Toy[iSample][iToy], samples[SampleObjectMap[iSample]]);
       chi2_pred_vec[iToy][iSample] = GetLLH(PredFluctHist, MC_Hist_Toy[iSample][iToy], W2_Hist_Toy[iSample][iToy], samples[SampleObjectMap[iSample]]);
  
       chi2_dat_vec[iToy].back()  += chi2_dat_vec[iToy][iSample];
       chi2_mc_vec[iToy].back()   += chi2_mc_vec[iToy][iSample];
       chi2_pred_vec[iToy].back() += chi2_pred_vec[iToy][iSample];
     }
   }
  
   MakeChi2Plots(chi2_mc_vec,   "-2LLH (Draw Fluc, Draw)", chi2_dat_vec, "-2LLH (Data, Draw)", SampleDir, "_drawfluc_draw");
   MakeChi2Plots(chi2_pred_vec, "-2LLH (Pred Fluc, Draw)", chi2_dat_vec, "-2LLH (Data, Draw)", SampleDir, "_predfluc_draw");
 }

◆ ProduceSpectra()

std::vector< std::unique_ptr< TH2D > > PredictiveThrower::ProduceSpectra	(	const std::vector< std::vector< std::unique_ptr< TH1D >>> &	Toys,
		const std::string	suffix
	)

private

Produce Violin style spectra.

Definition at line 438 of file PredictiveThrower.cpp.

                                                                                              {
 // *************************
   std::vector<double> MaxValue(TotalNumberOfSamples);
   #ifdef MULTITHREAD
   #pragma omp parallel for collapse(2)
   #endif
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     for (int toy = 0; toy < Ntoys; ++toy) {
       double max_val = Toys[sample][toy]->GetMaximum();
       MaxValue[sample] = std::max(MaxValue[sample], max_val);
     }
   }
  
   std::vector<std::unique_ptr<TH2D>> Spectra(TotalNumberOfSamples);
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     // Get MC histogram x-axis binning
     TH1D* refHist = Toys[sample][0].get();
  
     const int n_bins_x = refHist->GetNbinsX();
     std::vector<double> x_bin_edges(n_bins_x + 1);
     for (int b = 0; b <= n_bins_x; ++b) {
       x_bin_edges[b] = refHist->GetXaxis()->GetBinLowEdge(b + 1); // ROOT bins start at 1
     }
     // Last edge is upper edge of last bin:
     x_bin_edges[n_bins_x] = refHist->GetXaxis()->GetBinUpEdge(n_bins_x);
  
     constexpr int n_bins_y = 400;
     constexpr double y_min = 0.0;
     const double y_max = MaxValue[sample] * 1.05;
  
     // Create TH2D with variable binning on x axis
     Spectra[sample] = std::make_unique<TH2D>(
       (SampleNames[sample] + "_" + suffix).c_str(),   // name
       (SampleNames[sample] + "_" + suffix).c_str(),   // title
       n_bins_x, x_bin_edges.data(),                   // x axis bins
       n_bins_y, y_min, y_max                          // y axis bins
     );
  
     Spectra[sample]->SetDirectory(nullptr);
     Spectra[sample]->Sumw2(true);
   }
  
   #ifdef MULTITHREAD
   #pragma omp parallel for
   #endif
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     for (int toy = 0; toy < Ntoys; ++toy) {
       FastViolinFill(Spectra[sample].get(), Toys[sample][toy].get());
     }
   }
  
   return Spectra;
 }

◆ ProduceToys()

void PredictiveThrower::ProduceToys ( )

Produce toys by throwing from MCMC.

If we found toys then skip process of making new toys

Setup useful information for toy generation

this store value of parameters sampled from a chain

Definition at line 281 of file PredictiveThrower.cpp.

                                     {
 // *************************
   if(LoadToys()) return;
  
   SetupToyGeneration();
  
   auto PosteriorFileName = Get<std::string>(fitMan->raw()["Predictive"]["PosteriorFile"], __FILE__, __LINE__);
  
   MACH3LOG_INFO("Starting {}", __func__);
  
   outputFile->cd();
   double Penalty = 0, Weight = 1;
   int Draw = 0;
  
   TTree *ToyTree = new TTree("ToySummary", "ToySummary");
   ToyTree->Branch("Penalty", &Penalty, "Penalty/D");
   ToyTree->Branch("Weight", &Weight, "Weight/D");
   ToyTree->Branch("Draw", &Draw, "Draw/I");
   ToyTree->Branch("NModelParams", &NModelParams, "NModelParams/I");
  
   TDirectory* ToyDirectory = outputFile->mkdir("Toys");
   ToyDirectory->cd();
  
   for (size_t iPDF = 0; iPDF < samples.size(); iPDF++)
   {
     auto* MaCh3Sample = dynamic_cast<SampleHandlerFD*>(samples[iPDF]);
     for (int SampleIndex = 0; SampleIndex < MaCh3Sample->GetNsamples(); ++SampleIndex)
     {
       // Get nominal spectra and event rates
       TH1D* DataHist1D = static_cast<TH1D*>(MaCh3Sample->GetDataHist(1));
       Data_Hist[iPDF] = M3::Clone(DataHist1D, MaCh3Sample->GetTitle() + "_data");
       Data_Hist[iPDF]->Write((MaCh3Sample->GetTitle() + "_data").c_str());
  
       TH1D* MCHist1D = static_cast<TH1D*>(MaCh3Sample->GetMCHist(1));
       MC_Nom_Hist[iPDF] = M3::Clone(MCHist1D, MaCh3Sample->GetTitle() + "_mc");
       MCHist1D->Write((MaCh3Sample->GetTitle() + "_mc").c_str());
  
       TH1D* W2Hist1D = static_cast<TH1D*>(MaCh3Sample->GetW2Hist(1));
       W2_Nom_Hist[iPDF] = M3::Clone(W2Hist1D, MaCh3Sample->GetTitle() + "_w2");
       W2Hist1D->Write((MaCh3Sample->GetTitle() + "_w2").c_str());
       delete W2Hist1D;
     }
   }
  
   std::vector<std::vector<double>> branch_vals(systematics.size());
   std::vector<std::vector<std::string>> branch_name(systematics.size());
  
   TChain* PosteriorFile = new TChain("posteriors");
   PosteriorFile->Add(PosteriorFileName.c_str());
   unsigned int Step = 0;
   PosteriorFile->SetBranchAddress("step", &Step);
  
   if (PosteriorFile->GetBranch("Weight")) {
     PosteriorFile->SetBranchStatus("Weight", true);
     PosteriorFile->SetBranchAddress("Weight", &Weight);
   } else {
     MACH3LOG_WARN("Not applying reweighting weight");
     Weight = 1.0;
   }
  
   for (size_t s = 0; s < systematics.size(); ++s) {
     systematics[s]->MatchMaCh3OutputBranches(PosteriorFile, branch_vals[s], branch_name[s]);
   }
   //Get the burn-in from the config
   auto burn_in = Get<unsigned int>(fitMan->raw()["Predictive"]["BurnInSteps"], __FILE__, __LINE__);
  
   //DL: Adding sanity check for chains shorter than burn in
   const unsigned int maxNsteps = static_cast<unsigned int>(PosteriorFile->GetMaximum("step"));
   if(burn_in >= maxNsteps)
   {
     MACH3LOG_ERROR("You are running on a chain shorter than burn in cut");
     MACH3LOG_ERROR("Maximal value of nSteps: {}, burn in cut {}", maxNsteps, burn_in);
     MACH3LOG_ERROR("You will run into infinite loop");
     MACH3LOG_ERROR("You can make new chain or modify burn in cut");
     throw MaCh3Exception(__FILE__,__LINE__);
   }
  
   TStopwatch TempClock;
   TempClock.Start();
   for(int i = 0; i < Ntoys; i++)
   {
     if( i % (Ntoys/10) == 0) {
       MaCh3Utils::PrintProgressBar(i, Ntoys);
     }
     int entry = 0;
     Step = 0;
  
     //YSP: Ensures you get an entry from the mcmc even when burn_in is set to zero (Although not advised :p ).
     //Take 200k burn in steps, WP: Eb C in 1st peaky
     // If we have combined chains by hadd need to check the step in the chain
     // Note, entry is not necessarily same as step due to merged ROOT files, so can't choose entry in the range BurnIn - nEntries :(
     while(Step < burn_in){
       entry = random->Integer(static_cast<unsigned int>(PosteriorFile->GetEntries()));
       PosteriorFile->GetEntry(entry);
     }
     if(!Is_PriorPredictive) Draw = entry;
     for (size_t s = 0; s < systematics.size(); ++s)
     {
       systematics[s]->SetParameters(branch_vals[s]);
  
       //KS: Below line can help you get prior predictive distributions which are helpful for getting pre and post ND fit spectra
       //YSP: If not set in the config, the code runs SK Posterior Predictive distributions by default. If true, then the code runs SK prior predictive.
       if(Is_PriorPredictive) systematics[s]->ThrowParameters();
     }
  
     // This set some params to prior value this way you can evaluate errors from subset of errors
     SetParamters();
  
     Penalty = 0;
     if(FullLLH) {
       for (size_t s = 0; s < systematics.size(); ++s) {
         //KS: do times 2 because banff reports chi2
         Penalty = 2.0 * systematics[s]->GetLikelihood();
       }
     }
  
     PenaltyTerm[i] = Penalty;
     ReweightWeight[i] = Weight;
  
     for (size_t iPDF = 0; iPDF < samples.size(); iPDF++) {
       samples[iPDF]->Reweight();
     }
  
     for (size_t iPDF = 0; iPDF < samples.size(); iPDF++)
     {
       auto* MaCh3Sample = dynamic_cast<SampleHandlerFD*>(samples[iPDF]);
       for (int SampleIndex = 0; SampleIndex < MaCh3Sample->GetNsamples(); ++SampleIndex)
       {
         TH1D* MCHist1D = static_cast<TH1D*>(MaCh3Sample->GetMCHist(1));
         MC_Hist_Toy[iPDF][i] = M3::Clone(MCHist1D, MaCh3Sample->GetTitle() + "_mc_" + std::to_string(i));
         MC_Hist_Toy[iPDF][i]->Write();
  
         TH1D* W2Hist1D = static_cast<TH1D*>(MaCh3Sample->GetW2Hist(1));
         W2_Hist_Toy[iPDF][i] = M3::Clone(W2Hist1D, MaCh3Sample->GetTitle() + "_w2_" + std::to_string(i));
         W2_Hist_Toy[iPDF][i]->Write();
         delete W2Hist1D;
       }
     }
     ToyTree->Fill();
   }//end of toys loop
   TempClock.Stop();
  
   delete PosteriorFile;
   ToyDirectory->Close();
   delete ToyDirectory;
  
   outputFile->cd();
   ToyTree->Write();
   delete ToyTree;
  
   MACH3LOG_INFO("{} took {:.2f}s to finish for {} toys", __func__, TempClock.RealTime(), Ntoys);
 }

◆ RunMCMC()

void PredictiveThrower::RunMCMC ( )

inlineoverridevirtual

This is not used in this class.

Implements FitterBase.

Definition at line 38 of file PredictiveThrower.h.

                           {
     MACH3LOG_ERROR("{} is not supported in {}", __func__, GetName());
     throw MaCh3Exception(__FILE__ , __LINE__ );
   };

◆ RunPredictiveAnalysis()

void PredictiveThrower::RunPredictiveAnalysis ( )

Main routine responsible for producing posterior predictive distributions and $p$-value.

Definition at line 534 of file PredictiveThrower.cpp.

                                               {
 // *************************
   MACH3LOG_INFO("Starting {}", __func__);
   TStopwatch TempClock;
   TempClock.Start();
  
   auto DebugHistograms = GetFromManager<bool>(fitMan->raw()["Predictive"]["DebugHistograms"], false, __FILE__, __LINE__);
  
   TDirectory* PredictiveDir = outputFile->mkdir("Predictive");
   std::vector<TDirectory*> SampleDirectories;
   SampleDirectories.resize(TotalNumberOfSamples+1);
  
   for (int sample = 0; sample < TotalNumberOfSamples+1; ++sample) {
     SampleDirectories[sample] = PredictiveDir->mkdir(SampleNames[sample].c_str());
   }
  
   std::vector<std::unique_ptr<TH2D>> Spectra_mc = ProduceSpectra(MC_Hist_Toy, "mc");
   //std::vector<std::unique_ptr<TH2D>> Spectra_w2 = ProduceSpectra(W2_Hist_Toy, "w2");
   std::vector<std::unique_ptr<TH1D>> PostPred_mc(TotalNumberOfSamples);
   //std::vector<std::unique_ptr<TH1D>> PostPred_w2(TotalNumberOfSamples);
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     SampleDirectories[sample]->cd();
     Spectra_mc[sample]->Write();
     //Spectra_w2[sample]->Write();
  
     PostPred_mc[sample] = MakePredictive(MC_Hist_Toy[sample], SampleNames[sample], "mc", DebugHistograms);
     //PostPred_w2[sample] = MakePredictive(W2_Hist_Toy[sample], SampleNames[sample], "w2", DebugHistograms);
   }
  
   PosteriorPredictivepValue(PostPred_mc,
                             //PostPred_w2,
                             SampleDirectories);
  
   // Close directories
   for (int sample = 0; sample < TotalNumberOfSamples+1; ++sample) {
     SampleDirectories[sample]->Close();
     delete SampleDirectories[sample];
   }
  
   PredictiveDir->Close();
   delete PredictiveDir;
  
   outputFile->cd();
  
   TempClock.Stop();
   MACH3LOG_INFO("{} took {:.2f}s to finish for {} toys", __func__, TempClock.RealTime(), Ntoys);
 }

◆ SetParamters()

void PredictiveThrower::SetParamters ( )

private

This set some params to prior value this way you can evaluate errors from subset of errors.

Have ability to not throw legacy matrices

Alternatively vary only selected params

Definition at line 34 of file PredictiveThrower.cpp.

                                      {
 // *************************
   // WARNING This should be removed in the future
   auto DoNotThrowLegacyCov = GetFromManager<std::vector<std::string>>(fitMan->raw()["Predictive"]["DoNotThrowLegacyCov"], {}, __FILE__, __LINE__);
   for (size_t i = 0; i < DoNotThrowLegacyCov.size(); ++i) {
     for (size_t s = 0; s < systematics.size(); ++s) {
       if (systematics[s]->GetName() == DoNotThrowLegacyCov[i]) {
         systematics[s]->SetParameters();
         break;
       }
     }
   }
  
   if(ModelSystematic && ParameterGroupsNotVaried.size() > 0) ModelSystematic->SetGroupOnlyParameters(ParameterGroupsNotVaried);
  
   if (ModelSystematic && !ParameterOnlyToVary.empty()) {
     for (int i = 0; i < ModelSystematic->GetNumParams(); ++i) {
       if (ParameterOnlyToVary.find(i) == ParameterOnlyToVary.end()) {
         ModelSystematic->SetParProp(i, ModelSystematic->GetParInit(i));
       }
     }
   }
 }

◆ SetupSampleInformation()

void PredictiveThrower::SetupSampleInformation ( )

private

Setup sample information.

Definition at line 61 of file PredictiveThrower.cpp.

                                                {
 // *************************
   TotalNumberOfSamples = 0;
   for (size_t iPDF = 0; iPDF < samples.size(); iPDF++)
   {
     auto* MaCh3Sample = dynamic_cast<SampleHandlerFD*>(samples[iPDF]);
     if (!MaCh3Sample) {
       MACH3LOG_ERROR("Sample {} do not inherit from SampleHandlerFD this is not implemented", samples[iPDF]->GetTitle());
       throw MaCh3Exception(__FILE__, __LINE__);
     }
     TotalNumberOfSamples += samples[iPDF]->GetNsamples();
     if(samples[iPDF]->GetNsamples() > 1){
       MACH3LOG_ERROR("Sample has more than one sample {} ::", samples[iPDF]->GetNsamples());
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
   }
  
   MC_Hist_Toy.resize(TotalNumberOfSamples);
   W2_Hist_Toy.resize(TotalNumberOfSamples);
   Data_Hist.resize(TotalNumberOfSamples);
   MC_Nom_Hist.resize(TotalNumberOfSamples);
   W2_Nom_Hist.resize(TotalNumberOfSamples);
  
   SampleObjectMap.resize(TotalNumberOfSamples);
   SampleNames.resize(TotalNumberOfSamples+1);
  
   int currentIndex = 0;
   for (size_t iPDF = 0; iPDF < samples.size(); ++iPDF) {
     for (int subSampleIndex = 0; subSampleIndex < samples[iPDF]->GetNsamples(); ++subSampleIndex) {
       SampleObjectMap[currentIndex] = static_cast<int>(iPDF); // map the current global sample index to this sample object
       ++currentIndex;
     }
   }
  
   for (int sample = 0; sample < TotalNumberOfSamples; ++sample) {
     MC_Hist_Toy[sample].resize(Ntoys);
     W2_Hist_Toy[sample].resize(Ntoys);
     SampleNames[sample] = samples[sample]->GetTitle();
   }
   SampleNames[TotalNumberOfSamples] = "Total";
 }

◆ SetupToyGeneration()

void PredictiveThrower::SetupToyGeneration ( )

private

Setup useful variables etc before stating toy generation.

Let's ask the manager what are the file with covariance matrix

Definition at line 105 of file PredictiveThrower.cpp.

                                            {
 // *************************
   int counter = 0;
   for (size_t s = 0; s < systematics.size(); ++s) {
     auto* MaCh3Params = dynamic_cast<ParameterHandlerGeneric*>(systematics[s]);
     if(MaCh3Params) {
       ModelSystematic = MaCh3Params;
       counter++;
     }
   }
  
   SetupSampleInformation();
  
   if(Is_PriorPredictive) {
     MACH3LOG_INFO("You've chosen to run Prior Predictive Distribution");
   } else {
     auto PosteriorFileName = Get<std::string>(fitMan->raw()["Predictive"]["PosteriorFile"], __FILE__, __LINE__);
     //KS: We use MCMCProcessor to get names of covariances that were actually used to produce given chain
     MCMCProcessor Processor(PosteriorFileName);
     Processor.Initialise();
  
     YAML::Node ConfigInChain = Processor.GetCovConfig(kXSecPar);
     if(ModelSystematic){
       YAML::Node ConfigNow = ModelSystematic->GetConfig();
       if (!compareYAMLNodes(ConfigNow, ConfigInChain))
       {
         MACH3LOG_ERROR("Yaml configs in previous chain and current one are different", PosteriorFileName);
         throw MaCh3Exception(__FILE__ , __LINE__ );
       }
     }
   }
   if(counter > 1) {
     MACH3LOG_ERROR("Found {} ParmaterHandler inheriting from ParameterHandlerGeneric, I can accept at most 1", counter);
     throw MaCh3Exception(__FILE__, __LINE__);
   }
  
   for (size_t s = 0; s < systematics.size(); ++s) {
     NModelParams += systematics[s]->GetNumParams();
   }
  
   if (ModelSystematic) {
     auto ThrowParamGroupOnly = GetFromManager<std::vector<std::string>>(fitMan->raw()["Predictive"]["ThrowParamGroupOnly"], {}, __FILE__, __LINE__);
     auto UniqueParamGroup = ModelSystematic->GetUniqueParameterGroups();
     auto ParameterOnlyToVaryString = GetFromManager<std::vector<std::string>>(fitMan->raw()["Predictive"]["ThrowSinlgeParams"], {}, __FILE__, __LINE__);
  
     if (!ThrowParamGroupOnly.empty() && !ParameterOnlyToVaryString.empty()) {
       MACH3LOG_ERROR("Can't use ThrowParamGroupOnly and ThrowSinlgeParams at the same time");
       throw MaCh3Exception(__FILE__, __LINE__);
     }
  
     if (!ParameterOnlyToVaryString.empty()) {
       MACH3LOG_INFO("I will throw only: {}", fmt::join(ParameterOnlyToVaryString, ", "));
       std::vector<int> ParameterVary(ParameterOnlyToVaryString.size());
  
       for (size_t i = 0; i < ParameterOnlyToVaryString.size(); ++i) {
         ParameterVary[i] = ModelSystematic->GetParIndex(ParameterOnlyToVaryString[i]);
         if (ParameterVary[i] == M3::_BAD_INT_) {
           MACH3LOG_ERROR("Can't proceed if param {} is missing", ParameterOnlyToVaryString[i]);
           throw MaCh3Exception(__FILE__, __LINE__);
         }
       }
       ParameterOnlyToVary = std::unordered_set<int>(ParameterVary.begin(), ParameterVary.end());
     } else {
       MACH3LOG_INFO("I have following parameter groups: {}", fmt::join(UniqueParamGroup, ", "));
       if (ThrowParamGroupOnly.empty()) {
         MACH3LOG_INFO("I will vary all");
       } else {
         std::unordered_set<std::string> throwOnlySet(ThrowParamGroupOnly.begin(), ThrowParamGroupOnly.end());
         ParameterGroupsNotVaried.clear();
  
         for (const auto& group : UniqueParamGroup) {
           if (throwOnlySet.find(group) == throwOnlySet.end()) {
             ParameterGroupsNotVaried.push_back(group);
           }
         }
  
         MACH3LOG_INFO("I will vary: {}", fmt::join(ThrowParamGroupOnly, ", "));
         MACH3LOG_INFO("Exclude: {}", fmt::join(ParameterGroupsNotVaried, ", "));
       }
     }
   }
 }

Member Data Documentation

◆ Data_Hist

std::vector<std::unique_ptr<TH1D> > PredictiveThrower::Data_Hist

private

Vector of Data histograms.

Definition at line 114 of file PredictiveThrower.h.

◆ FullLLH

bool PredictiveThrower::FullLLH

private

KS: Use Full LLH or only sample contribution based on discussion with Asher we almost always only want the sample likelihood.

Definition at line 90 of file PredictiveThrower.h.

◆ Is_PriorPredictive

bool PredictiveThrower::Is_PriorPredictive

private

Whether it is Prior or Posterior predictive.

Definition at line 94 of file PredictiveThrower.h.

◆ MC_Hist_Toy

std::vector<std::vector<std::unique_ptr<TH1D> > > PredictiveThrower::MC_Hist_Toy

private

Vector of MC histograms per sample and toy experiment. Indexed as [sample][toy].

Definition at line 122 of file PredictiveThrower.h.

◆ MC_Nom_Hist

std::vector<std::unique_ptr<TH1D> > PredictiveThrower::MC_Nom_Hist

private

Vector of MC histograms.

Definition at line 116 of file PredictiveThrower.h.

◆ ModelSystematic

ParameterHandlerGeneric* PredictiveThrower::ModelSystematic

private

Pointer to El Generico.

Definition at line 111 of file PredictiveThrower.h.

◆ NModelParams

int PredictiveThrower::NModelParams

private

KS: Count total number of model parameters which can be used for stuff like BIC.

Definition at line 92 of file PredictiveThrower.h.

◆ Ntoys

int PredictiveThrower::Ntoys

private

Number of toys we are generating analysing.

Definition at line 104 of file PredictiveThrower.h.

◆ ParameterGroupsNotVaried

std::vector<std::string> PredictiveThrower::ParameterGroupsNotVaried

private

KS: Names of parameter groups that will not be varied.

Definition at line 106 of file PredictiveThrower.h.

◆ ParameterOnlyToVary

std::unordered_set<int> PredictiveThrower::ParameterOnlyToVary

private

KS: Index of parameters groups that will be varied.

Definition at line 108 of file PredictiveThrower.h.

◆ PenaltyTerm

std::vector<double> PredictiveThrower::PenaltyTerm

private

Penalty term values for each toy by default 0.

Definition at line 130 of file PredictiveThrower.h.

◆ ReweightWeight

std::vector<double> PredictiveThrower::ReweightWeight

private

Reweighting factors applied for each toy, by default 1.

Definition at line 128 of file PredictiveThrower.h.

◆ SampleNames

std::vector<std::string> PredictiveThrower::SampleNames

private

Name of a single sample.

Definition at line 99 of file PredictiveThrower.h.

◆ SampleObjectMap

std::vector<int> PredictiveThrower::SampleObjectMap

private

Maps if sample with given SampleHandler, useful if we have more than one sample in single object.

Definition at line 101 of file PredictiveThrower.h.

◆ TotalNumberOfSamples

int PredictiveThrower::TotalNumberOfSamples

private

Number of toys we are generating analysing.

Definition at line 97 of file PredictiveThrower.h.

◆ W2_Hist_Toy

std::vector<std::vector<std::unique_ptr<TH1D> > > PredictiveThrower::W2_Hist_Toy

private

Vector of W² histograms per sample and toy experiment. Indexed as [sample][toy]

Definition at line 125 of file PredictiveThrower.h.

◆ W2_Nom_Hist

std::vector<std::unique_ptr<TH1D> > PredictiveThrower::W2_Nom_Hist

private

Vector of W2 histograms.

Definition at line 118 of file PredictiveThrower.h.

The documentation for this class was generated from the following files:

Fitters/PredictiveThrower.h
Fitters/PredictiveThrower.cpp

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ PredictiveThrower()

◆ ~PredictiveThrower()

Member Function Documentation

◆ GetLLH()

◆ LoadToys()

◆ MakeChi2Plots()

◆ MakePredictive()

◆ PosteriorPredictivepValue()

◆ ProduceSpectra()

◆ ProduceToys()

◆ RunMCMC()

◆ RunPredictiveAnalysis()

◆ SetParamters()

◆ SetupSampleInformation()

◆ SetupToyGeneration()

Member Data Documentation

◆ Data_Hist

◆ FullLLH

◆ Is_PriorPredictive

◆ MC_Hist_Toy

◆ MC_Nom_Hist

◆ ModelSystematic

◆ NModelParams

◆ Ntoys

◆ ParameterGroupsNotVaried

◆ ParameterOnlyToVary

◆ PenaltyTerm

◆ ReweightWeight

◆ SampleNames

◆ SampleObjectMap

◆ TotalNumberOfSamples

◆ W2_Hist_Toy

◆ W2_Nom_Hist