MonteCarloEngine Class Reference

An engine class for Monte Carlo. More...

#include <MonteCarloEngine.h>

Inheritance diagram for MonteCarloEngine:
[legend]
Collaboration diagram for MonteCarloEngine:
[legend]

Detailed Description

An engine class for Monte Carlo.

Author
HEPfit Collaboration

The MonteCarlo engine is has overloaded from BCEngineMCMC class in the BAT libraries and user defined fumctions to facilitate the Monte Carlo run.

Definition at line 46 of file MonteCarloEngine.h.

Public Member Functions

void AddChains ()
 A method to add the observable values and weights to the chain information. More...
 
void CheckHistogram (const TH1D &hist, const std::string name)
 This member checks if there is overflow of the 1D histogram. More...
 
void CheckHistogram (const TH2D &hist, const std::string name)
 This member checks if there is overflow of the 2D histogram. More...
 
double computeNormalization ()
 
std::string computeStatistics ()
 A get method to compute the mean and rms of the computed observables. More...
 
void DefineParameters ()
 A member to classify the prior of the model parameters varied in the Monte Carlo. More...
 
double FirstDerivative (const BCParameter *par, std::vector< double > point)
 
double Function_h (std::vector< double > point)
 
std::map< std::string, BCH1D * > getHistograms1D () const
 
std::map< std::string, BCH2D * > getHistograms2D () const
 
std::string getHistoLog () const
 A get method to access the stream that stores the log messages coming from histogram printing and checking. More...
 
int getNumOfDiscardedEvents () const
 A get method to access the number of events discarded due to failure to update model. These events are not used for the MCMC run. More...
 
int getNumOfUsedEvents () const
 A get method to access the number of events used in the MCMC run. More...
 
void Initialize (Model *Mod_i)
 Initialization of the Monte Carlo Engine. More...
 
double LogLikelihood (const std::vector< double > &parameters)
 This member calculates the loglikelihood for the observables included in the MCMC run. More...
 
void MCMCIterationInterface ()
 Overloaded from BCEngineMCMC in BAT More...
 
 MonteCarloEngine (const std::vector< ModelParameter > &ModPars_i, boost::ptr_vector< Observable > &Obs_i, std::vector< Observable2D > &Obs2D_i, std::vector< CorrelatedGaussianObservables > &CGO_i, std::vector< CorrelatedGaussianParameters > &CGP_i)
 Constructor. More...
 
void PrintCorrelationMatrix (const std::string filename)
 This member generates the correlation matrix using BCH2D from the BAT libraries. More...
 
void PrintHistogram (BCModelOutput &out, Observable &it, const std::string OutputDir)
 Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables. More...
 
void PrintHistogram (BCModelOutput &out, const std::string OutputDir)
 Member used for printing histograms for observables, observable2D, correlated Gaussian observables and model parameters vs. observables. More...
 
double SecondDerivative (const BCParameter *par1, const BCParameter *par2, std::vector< double > point)
 
void setDParsFromParameters (const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
 
void setHistogramOverFlow (bool overflow)
 
void setNChains (unsigned int i)
 A set method to fix the number of chains. More...
 
std::string writePreRunData ()
 A method to write in a text file the best fit parameters and the prerun scale factors. More...
 
 ~MonteCarloEngine ()
 The default destructor. Some pointers defined in this class are explicitly freed. More...
 

Private Attributes

std::vector< CorrelatedGaussianObservables > & CGO
 A vector of correlated Gaussian observables. More...
 
const std::vector< CorrelatedGaussianParameters > & CGP
 A vector of correlated Gaussian parameters. More...
 
std::map< std::string, TPrincipal * > CorrelationMap
 A map between the name of a theory observable and its maximum value. More...
 
std::map< std::string, double > DPars
 A map between parameter names and their values. More...
 
std::map< std::string, BCH1D * > Histo1D
 A map between pointers to objects of type BCH1D (BAT) and their given names. More...
 
std::map< std::string, BCH2D * > Histo2D
 A map between pointers to objects of type BCH2D (BAT) and their given names. More...
 
std::ostringstream HistoLog
 A stream to store the output messages from printing and checking histograms. More...
 
unsigned int kmax
 The number of observables. More...
 
unsigned int kwmax
 The number of observables whose weights are used for the MCMC. More...
 
ModelMod
 A pointer to an abject of type Model. More...
 
const std::vector< ModelParameter > & ModPars
 A vector of model parameters. More...
 
int NumOfDiscardedEvents
 The number of events for which the update of the model fails and these events are not used for the MCMC run. More...
 
int NumOfUsedEvents
 The number of events for which the model is successfully updated and hence used for the MCMC run. More...
 
std::vector< Observable2D > & Obs2D_ALL
 A vector of all pairs of observable for Observable2D. More...
 
boost::ptr_vector< Observable > & Obs_ALL
 A vector of all observables. More...
 
double * obval
 A pointer to the vector of observable values. More...
 
double * obweight
 A pointer to the vector of observable weights. More...
 
int rank
 Rank of the process for a MPI run. Value is 0 for a serial run. More...
 
std::map< std::string, double > thMax
 A map between the name of a theory observable and its maximum value. More...
 
std::map< std::string, double > thMin
 A map between the name of a theory observable and its minimum value. More...
 

Constructor & Destructor Documentation

MonteCarloEngine::MonteCarloEngine ( const std::vector< ModelParameter > &  ModPars_i,
boost::ptr_vector< Observable > &  Obs_i,
std::vector< Observable2D > &  Obs2D_i,
std::vector< CorrelatedGaussianObservables > &  CGO_i,
std::vector< CorrelatedGaussianParameters > &  CGP_i 
)

Constructor.

Parameters
[in]ModPars_ithe vector of model parameters as defined in SomeModel.conf
[in]Obs_ithe vector of observables as defined in SomeModel.conf
[in]Obs2D_ithe vector of observables2D as defined in SomeModel.conf
[in]CGO_ithe vector of correlated Gaussian observables as defined in SomeModel.conf

Definition at line 23 of file MonteCarloEngine.cpp.

29 : BCModel(""), ModPars(ModPars_i), CGP(CGP_i), Obs_ALL(Obs_i), Obs2D_ALL(Obs2D_i),
31  obval = NULL;
32  obweight = NULL;
33  Mod = NULL;
34  TH1::StatOverflows(kTRUE);
35 #ifdef _MPI
36  rank = MPI::COMM_WORLD.Get_rank();
37 #else
38  rank = 0;
39 #endif
40 };
int NumOfDiscardedEvents
The number of events for which the update of the model fails and these events are not used for the MC...
const std::vector< ModelParameter > & ModPars
A vector of model parameters.
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
const std::vector< CorrelatedGaussianParameters > & CGP
A vector of correlated Gaussian parameters.
double * obval
A pointer to the vector of observable values.
int NumOfUsedEvents
The number of events for which the model is successfully updated and hence used for the MCMC run...
int rank
Rank of the process for a MPI run. Value is 0 for a serial run.
double * obweight
A pointer to the vector of observable weights.
Model * Mod
A pointer to an abject of type Model.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
MonteCarloEngine::~MonteCarloEngine ( )

The default destructor. Some pointers defined in this class are explicitly freed.

Definition at line 132 of file MonteCarloEngine.cpp.

134 {
135  delete [] obval;
136  delete [] obweight;
137  /* The following code has been commented out pending further review.
138  It is causing crashes at the termination of the code if the histograms
139  are accessed from the main program.*/
140  // for (std::map<std::string, BCH1D *>::iterator it = Histo1D.begin();
141  // it != Histo1D.end(); it++)
142  // delete it->second;
143  // for (std::map<std::string, BCH2D *>::iterator it = Histo2D.begin();
144  // it != Histo2D.end(); it++)
145  // delete it->second;
146  for (std::map<std::string, TPrincipal *>::iterator it = CorrelationMap.begin();
147  it != CorrelationMap.end(); it++)
148  delete it->second;
149 };
double * obval
A pointer to the vector of observable values.
std::map< std::string, TPrincipal * > CorrelationMap
A map between the name of a theory observable and its maximum value.
double * obweight
A pointer to the vector of observable weights.

Member Function Documentation

void MonteCarloEngine::AddChains ( )

A method to add the observable values and weights to the chain information.

If the WriteChain flag is set to true in the MonteCarlo.conf file then chain information is written. This method adds the observable values and weights (for the observables included in the MCMC) to the chain information.

Definition at line 590 of file MonteCarloEngine.cpp.

590  {
591  fMCMCFlagWritePreRunToFile = false;
592  int k = 0, kweight = 0;
593  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
594  it < Obs_ALL.end(); it++) {
595  if (!it->isTMCMC()) {
596  for (unsigned int i = 0; i < fMCMCNChains; ++i)
597  fMCMCTrees[i]->Branch(it->getName().c_str(), &obval[i * kmax + k],
598  (it->getName() + "/D").c_str());
599  k++;
600  if (it->getDistr().compare("noweight") != 0) {
601  for (unsigned int i = 0; i < fMCMCNChains; ++i)
602  fMCMCTrees[i]->Branch((it->getName() + "_weight").c_str(),
603  &obweight[i * kwmax + kweight],
604  (it->getName() + "_weight/D").c_str());
605  kweight++;
606  }
607  }
608  }
609 }
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
unsigned int kwmax
The number of observables whose weights are used for the MCMC.
double * obval
A pointer to the vector of observable values.
double * obweight
A pointer to the vector of observable weights.
unsigned int kmax
The number of observables.
void MonteCarloEngine::CheckHistogram ( const TH1D &  hist,
const std::string  name 
)

This member checks if there is overflow of the 1D histogram.

Parameters
[in]hista reference to an object of type TH1D as defined in the ROOT libraries
[in]namethe name for the histogram

Definition at line 494 of file MonteCarloEngine.cpp.

494  {
495  // output the portions of underflow and overflow bins
496  double UnderFlowContent = hist.GetBinContent(0);
497  double OverFlowContent = hist.GetBinContent(NBINS1D + 1);
498  double Integral = hist.Integral();
499  double TotalContent = 0.0;
500  for (int n = 0; n <= NBINS1D + 1; n++)
501  TotalContent += hist.GetBinContent(n);
502  HistoLog << name << ": "
503  << Integral / TotalContent * 100. << "% within the range, "
504  << UnderFlowContent / TotalContent * 100. << "% underflow, "
505  << OverFlowContent / TotalContent * 100. << "% overflow"
506  << std::endl;
507 }
#define NBINS1D
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
void MonteCarloEngine::CheckHistogram ( const TH2D &  hist,
const std::string  name 
)

This member checks if there is overflow of the 2D histogram.

Parameters
[in]hista reference to an object of type TH2D as defined in the ROOT libraries
[in]namethe name for the histogram

Definition at line 509 of file MonteCarloEngine.cpp.

509  {
510  double Integral = hist.Integral();
511  double TotalContent = 0.0;
512  for (int m = 0; m <= NBINS2D + 1; m++)
513  for (int n = 0; n <= NBINS2D + 1; n++)
514  TotalContent += hist.GetBinContent(m, n);
515  HistoLog << name << ": "
516  << Integral / TotalContent * 100. << "% within the ranges"
517  << std::endl;
518 }
#define NBINS2D
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
double MonteCarloEngine::computeNormalization ( )

Definition at line 832 of file MonteCarloEngine.cpp.

832  {
833 
834  unsigned int Npars = GetNParameters();
835  std::vector<double> mode(GetBestFitParameters());
836  gslpp::matrix<double> Hessian(Npars, Npars, 0.);
837 
838  for (unsigned int i = 0; i < Npars; i++)
839  for (unsigned int j = 0; j < Npars; j++) {
840  // calculate Hessian matrix element
841  Hessian.assign(i, j, -SecondDerivative(GetParameter(i), GetParameter(j), GetBestFitParameters()));
842  }
843 
844  double det_Hessian = Hessian.determinant();
845 
846  return exp(Npars / 2. * log(2. * M_PI) + 0.5 * log(1. / det_Hessian) + LogLikelihood(mode) + LogAPrioriProbability(mode));
847 }
A class for constructing and defining operations on real matrices.
double SecondDerivative(const BCParameter *par1, const BCParameter *par2, std::vector< double > point)
double LogLikelihood(const std::vector< double > &parameters)
This member calculates the loglikelihood for the observables included in the MCMC run...
complex log(const complex &z)
complex exp(const complex &z)
std::string MonteCarloEngine::computeStatistics ( )

A get method to compute the mean and rms of the computed observables.

Returns
a string containing the statistics

Definition at line 640 of file MonteCarloEngine.cpp.

640  {
641  std::vector<double> mode(GetBestFitParameters());
642  if (mode.size() == 0) {
643  if(rank == 0) throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. Observables statistics cannot be generated.\n");
644  else sleep (2);
645  }
646  std::ostringstream StatsLog;
647  int i = 0;
648  StatsLog << "Statistics file for Observables, Binned Observables and Corellated Gaussian Observables.\n" << std::endl;
649  if (Obs_ALL.size() > 0) StatsLog << "Observables:\n" << std::endl;
650  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
651 
652  if (it->getObsType().compare("BinnedObservable") == 0) {
653  StatsLog << " (" << ++i << ") Binned Observable \"";
654  StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << ", " << it->getTho()->getBinMax() << "]" << "\":";
655  } else if (it->getObsType().compare("FunctionObservable") == 0) {
656  StatsLog << " (" << ++i << ") Function Observable \"";
657  StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << "]" << "\":";
658  } else if (it->getObsType().compare("HiggsObservable") == 0) {
659  StatsLog << " (" << ++i << ") Higgs Observable \"";
660  StatsLog << it->getName() << "\":";
661  } else {
662  StatsLog << " (" << ++i << ") " << it->getObsType() << " \"";
663  StatsLog << it->getName() << "\":";
664  }
665 
666  StatsLog << std::endl;
667 
668  BCH1D * bch1d = Histo1D[it->getName()];
669  if (bch1d->GetHistogram()->Integral() > 0.0) {
670  StatsLog << " Mean +- sqrt(V): " << std::setprecision(4)
671  << bch1d->GetMean() << " +- " << std::setprecision(4)
672  << bch1d->GetRMS() << std::endl
673 
674  << " (Marginalized) mode: " << bch1d->GetMode() << std::endl;
675 
676  std::vector<double> v1;
677  std::vector<double> v2;
678  std::vector<double> v3;
679  v1 = bch1d->GetSmallestIntervals(0.6827);
680  v2 = bch1d->GetSmallestIntervals(0.9545);
681  v3 = bch1d->GetSmallestIntervals(0.9973);
682  StatsLog << " Smallest interval(s) containing at least 68.27% and local mode(s):"
683  << std::endl;
684  for (unsigned j = 0; j < v1.size(); j += 5)
685  StatsLog << " (" << v1[j] << ", " << v1[j + 1]
686  << ") (local mode at " << v1[j + 3] << " with rel. height "
687  << v1[j + 2] << "; rel. area " << v1[j + 4] << ")"
688  << std::endl;
689  StatsLog << std::endl;
690 
691  StatsLog << " Smallest interval(s) containing at least 95.45% and local mode(s):"
692  << std::endl;
693  for (unsigned j = 0; j < v2.size(); j += 5)
694  StatsLog << " (" << v2[j] << ", " << v2[j + 1]
695  << ") (local mode at " << v2[j + 3] << " with rel. height "
696  << v2[j + 2] << "; rel. area " << v2[j + 4] << ")"
697  << std::endl;
698  StatsLog << std::endl;
699 
700  StatsLog << " Smallest interval(s) containing at least 99.73% and local mode(s):"
701  << std::endl;
702  for (unsigned j = 0; j < v3.size(); j += 5)
703  StatsLog << " (" << v3[j] << ", " << v3[j + 1]
704  << ") (local mode at " << v3[j + 3] << " with rel. height "
705  << v3[j + 2] << "; rel. area " << v3[j + 4] << ")"
706  << std::endl;
707  StatsLog << std::endl;
708  } else {
709  StatsLog << "\nWARNING: The histogram of " << it->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
710  }
711  }
712 
713  if (CGO.size() > 0) StatsLog << "\nCorellated (Gaussian) Observables:\n" << std::endl;
714  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
715  it1 < CGO.end(); it1++) {
716  StatsLog << "\n" << it1->getName() << ":\n" << std::endl;
717  i = 0;
718  std::vector<Observable> CGObs(it1->getObs());
719  for (std::vector<Observable>::iterator it2 = CGObs.begin(); it2 < CGObs.end(); it2++) {
720 
721  if (it2->getObsType().compare("BinnedObservable") == 0) {
722  StatsLog << " (" << ++i << ") Binned Observable \"";
723  StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << ", " << it2->getTho()->getBinMax() << "]" << "\":";
724  } else if (it2->getObsType().compare("FunctionObservable") == 0) {
725  StatsLog << " (" << ++i << ") Function Observable \"";
726  StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << "]" << "\":";
727  } else if (it2->getObsType().compare("HiggsObservable") == 0) {
728  StatsLog << " (" << ++i << ") Higgs Observable \"";
729  StatsLog << it2->getName() << "\":";
730  } else {
731  StatsLog << " (" << ++i << ") " << it2->getObsType() << " \"";
732  StatsLog << it2->getName() << "\":";
733  }
734 
735  StatsLog << std::endl;
736  BCH1D * bch1d = Histo1D[it2->getName()];
737  if (bch1d->GetHistogram()->Integral() > 0.0) {
738  StatsLog << " Mean +- sqrt(V): " << std::setprecision(4)
739  << bch1d->GetMean() << " +- " << std::setprecision(4)
740  << bch1d->GetRMS() << std::endl
741 
742  << " (Marginalized) mode: " << bch1d->GetMode() << std::endl;
743 
744  std::vector<double> v1;
745  std::vector<double> v2;
746  std::vector<double> v3;
747  v1 = bch1d->GetSmallestIntervals(0.6827);
748  v2 = bch1d->GetSmallestIntervals(0.9545);
749  v3 = bch1d->GetSmallestIntervals(0.9973);
750  StatsLog << " Smallest interval(s) containing at least 68.27% and local mode(s):"
751  << std::endl;
752  for (unsigned j = 0; j < v1.size(); j += 5)
753  StatsLog << " (" << v1[j] << ", " << v1[j + 1]
754  << ") (local mode at " << v1[j + 3] << " with rel. height "
755  << v1[j + 2] << "; rel. area " << v1[j + 4] << ")"
756  << std::endl;
757  StatsLog << std::endl;
758 
759  StatsLog << " Smallest interval(s) containing at least 95.45% and local mode(s):"
760  << std::endl;
761  for (unsigned j = 0; j < v2.size(); j += 5)
762  StatsLog << " (" << v2[j] << ", " << v2[j + 1]
763  << ") (local mode at " << v2[j + 3] << " with rel. height "
764  << v2[j + 2] << "; rel. area " << v2[j + 4] << ")"
765  << std::endl;
766  StatsLog << std::endl;
767 
768  StatsLog << " Smallest interval(s) containing at least 99.73% and local mode(s):"
769  << std::endl;
770  for (unsigned j = 0; j < v3.size(); j += 5)
771  StatsLog << " (" << v3[j] << ", " << v3[j + 1]
772  << ") (local mode at " << v3[j + 3] << " with rel. height "
773  << v3[j + 2] << "; rel. area " << v3[j + 4] << ")"
774  << std::endl;
775  StatsLog << std::endl;
776  } else {
777  StatsLog << "\nWARNING: The histogram of " << it2->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
778  }
779  }
780  if (it1->isPrediction()) {
781  int size = it1->getObs().size();
782  CorrelationMap[it1->getName()]->MakePrincipals();
783  //CorrelationMap[it1->getName()]->Print();
784  TMatrixD * corr = const_cast<TMatrixD*>(CorrelationMap[it1->getName()]->GetCovarianceMatrix());
785  *corr *= (double)size;
786  StatsLog << "\nThe correlation matrix for " << it1->getName() << " is given by the " << size << "x"<< size << " matrix:\n" << std::endl;
787 
788  for (int i = 0; i < size + 1; i++) {
789  if (i == 0) StatsLog << std::setw(4) << "" << " | ";
790  else StatsLog << std::setw(5) << i << std::setw(5) << " |";
791  }
792  StatsLog << std::endl;
793  for (int i = 0; i < size + 1; i++) {
794  if (i == 0) StatsLog << std::setw(8) << "--------";
795  else StatsLog << std::setw(10) << "----------";
796  }
797  StatsLog << std::endl;
798  for (int i = 0; i < size; i++) {
799  for (int j = 0; j < size + 1; j++) {
800  if (j == 0) StatsLog << std::setw(4) << i+1 << " |";
801  else StatsLog << std::setprecision(5) << std::setw(10) << (*corr)(i, j - 1);
802  }
803  StatsLog << std::endl;
804  }
805  }
806  }
807 
808  std::vector<double> parsa;
809  for (unsigned int i = 0; i < GetNParameters(); i++)
810  parsa.push_back(GetMarginalized(fParameters[i])->GetMean());
811  double llik = LogLikelihood(parsa);
812  StatsLog << "LogLikelihood on the mean values of parameters: " << llik << std::endl;
813  double llika = Histo1D["LogLikelihood"]->GetMean();
814  StatsLog << "LogLikelihood mean value: " << llika << std::endl;
815  double dbar = -2.*llika; //Wikipedia notation...
816  double dothetabar = -2.*llik; //Wikipedia notation...
817  StatsLog << "IC value (don't ask me what it means...): " << 3.*dbar - 2.*dothetabar << std::endl;
818  StatsLog << "DIC value (same as above...): " << 2.*dbar - dothetabar << std::endl;
819  return StatsLog.str().c_str();
820 }
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
double LogLikelihood(const std::vector< double > &parameters)
This member calculates the loglikelihood for the observables included in the MCMC run...
int rank
Rank of the process for a MPI run. Value is 0 for a serial run.
std::map< std::string, TPrincipal * > CorrelationMap
A map between the name of a theory observable and its maximum value.
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
void MonteCarloEngine::DefineParameters ( )

A member to classify the prior of the model parameters varied in the Monte Carlo.

The model parameters being varied are first sorted out checking for the existence of of Gaussian error ( \(\delta_{g}\)) or flat error ( \(\delta_{f}\)). The SetPriorGaus method and SetPriorConstant method in BCEngineMCMC (BAT) is used to set a Gaussian or flat prior for the model parameters respectively. In case a model parameter has both a Gaussian and a flat error then the combined function is built using TF1. The definition of the combined error is

\[ \delta_{combined}(x) = \frac{1}{4\delta_f}\left(\frac{Erf(x - x_{ave}+\delta_f)}{\sqrt{2}\delta_g} - \frac{Erf(x - x_{ave}-\delta_f)}{\sqrt{2}\delta_g}\right), \]

with

\[ Erf(x)=\frac{1}{\sqrt{2}\pi}\int_0^x \exp^{-\frac{t^2}{2}}dt. \]

Definition at line 153 of file MonteCarloEngine.cpp.

153  {
154  // Add parameters to your model here.
155  // You can then use them in the methods below by calling the
156  // parameters.at(i) or parameters[i], where i is the index
157  // of the parameter. The indices increase from 0 according to the
158  // order of adding the parameters.
159  if (rank == 0) std::cout << "\nParameters varied in Monte Carlo:" << std::endl;
160  int k = 0;
161  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin();
162  it < ModPars.end(); it++) {
163  if (it->geterrf() == 0. && it->geterrg() == 0.)
164  continue;
165  AddParameter(it->getname().c_str(), it->getmin(), it->getmax());
166  if (rank == 0) std::cout << k << ": " << it->getname() << ", ";
167  if (it->IsCorrelated()) {
168  for (unsigned int i = 0; i < CGP.size(); i++) {
169  if (CGP[i].getName().compare(it->getCgp_name()) == 0) {
170  std::string index = it->getname().substr(CGP[i].getName().size());
171  long int lindex = strtol(index.c_str(),NULL,10);
172  if (lindex > 0)
173  DPars[CGP[i].getPar(lindex - 1).getname()] = 0.;
174  else {
175  std::stringstream out;
176  out << it->getname();
177  throw std::runtime_error("MonteCarloEngine::DefineParameters(): " + out.str() + "seems to be part of a CorrelatedGaussianParameters object, but I couldn't find the corresponding object");
178  }
179  }
180  }
181  } else
182  DPars[it->getname()] = 0.;
183  if (it->geterrf() == 0.) SetPriorGauss(k, it->getave(), it->geterrg());
184  else if (it->geterrg() == 0.) SetPriorConstant(k);
185  else {
186  TF1 * combined = new TF1(it->getname().c_str(),
187  "(TMath::Erf((x-[0]+[2])/sqrt(2.)/[1])-TMath::Erf((x-[0]-[2])/sqrt(2.)/[1]))/4./[2]",
188  it->getmin(), it->getmax());
189  combined->SetParameter(0, it->getave());
190  combined->SetParameter(1, it->geterrg());
191  combined->SetParameter(2, it->geterrf());
192  SetPrior(k, combined);
193  delete combined;
194  }
195  k++;
196  }
197 }
const std::vector< ModelParameter > & ModPars
A vector of model parameters.
const std::vector< CorrelatedGaussianParameters > & CGP
A vector of correlated Gaussian parameters.
std::map< std::string, double > DPars
A map between parameter names and their values.
int rank
Rank of the process for a MPI run. Value is 0 for a serial run.
double MonteCarloEngine::FirstDerivative ( const BCParameter *  par,
std::vector< double >  point 
)

Definition at line 884 of file MonteCarloEngine.cpp.

884  {
885 
886  if (point.size() != GetNParameters()) {
887  throw std::runtime_error("MCMCENgine::FirstDerivative : Invalid number of entries in the vector.");
888  }
889 
890  // define steps
891  const double dx1 = par->GetRangeWidth() / NSTEPS;
892  const double dx2 = dx1 * 2.;
893  const double dx3 = dx1 * 3.;
894 
895  // define points at which to evaluate
896  std::vector<double> x1p = point;
897  std::vector<double> x1m = point;
898  std::vector<double> x2p = point;
899  std::vector<double> x2m = point;
900  std::vector<double> x3p = point;
901  std::vector<double> x3m = point;
902 
903  unsigned idx = GetParameters().Index(par->GetName());
904 
905  x1p[idx] += dx1;
906  x1m[idx] -= dx1;
907  x2p[idx] += dx2;
908  x2m[idx] -= dx2;
909  x3p[idx] += dx3;
910  x3m[idx] -= dx3;
911 
912  const double m1 = (Function_h(x1p) - Function_h(x1m)) / 2. / dx1;
913  const double m2 = (Function_h(x2p) - Function_h(x2m)) / 4. / dx1;
914  const double m3 = (Function_h(x3p) - Function_h(x3m)) / 6. / dx1;
915 
916  return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
917 }
#define NSTEPS
double Function_h(std::vector< double > point)
double MonteCarloEngine::Function_h ( std::vector< double >  point)

Definition at line 919 of file MonteCarloEngine.cpp.

919  {
920  if (point.size() != GetNParameters()) {
921  throw std::runtime_error("MCMCENgine::Function_h : Invalid number of entries in the vector.");
922  }
923  return LogLikelihood(point) + LogAPrioriProbability(point);
924 }
double LogLikelihood(const std::vector< double > &parameters)
This member calculates the loglikelihood for the observables included in the MCMC run...
std::map<std::string, BCH1D * > MonteCarloEngine::getHistograms1D ( ) const
inline

Definition at line 225 of file MonteCarloEngine.h.

226  {
227  return Histo1D;
228  }
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
std::map<std::string, BCH2D * > MonteCarloEngine::getHistograms2D ( ) const
inline

Definition at line 230 of file MonteCarloEngine.h.

231  {
232  return Histo2D;
233  }
std::map< std::string, BCH2D * > Histo2D
A map between pointers to objects of type BCH2D (BAT) and their given names.
std::string MonteCarloEngine::getHistoLog ( ) const
inline

A get method to access the stream that stores the log messages coming from histogram printing and checking.

Returns
a string containing the log messages

Definition at line 182 of file MonteCarloEngine.h.

183  {
184  return HistoLog.str().c_str();
185  }
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
int MonteCarloEngine::getNumOfDiscardedEvents ( ) const
inline

A get method to access the number of events discarded due to failure to update model. These events are not used for the MCMC run.

Returns
the number of discarded events

Definition at line 211 of file MonteCarloEngine.h.

212  {
213  return NumOfDiscardedEvents;
214  }
int NumOfDiscardedEvents
The number of events for which the update of the model fails and these events are not used for the MC...
int MonteCarloEngine::getNumOfUsedEvents ( ) const
inline

A get method to access the number of events used in the MCMC run.

Returns
the number of events used in the MCMC run

Definition at line 220 of file MonteCarloEngine.h.

221  {
222  return NumOfUsedEvents;
223  }
int NumOfUsedEvents
The number of events for which the model is successfully updated and hence used for the MCMC run...
void MonteCarloEngine::Initialize ( Model Mod_i)

Initialization of the Monte Carlo Engine.

The initialization of the Monte Carlo Engine performs the following tasks

  • The observables and observables2D lists are checked and sorted into the ones which are included in the MCMC and the ones which are not.
  • The experimental likelihood for observables and observalbes2D are read if provided through a .root file and the corresponding output histograms are initiated.
  • The correlated Gaussian observables are sorted for inclusion in the MCMC and the corresponding output histograms are initiated.
  • The histograms for the model parameters vs. observables are initiated.
    The output specifications are passed from the SomeModel.conf file parsed by the InputParser. The number of bins for the 1D and 2D histograms are defined through the NBINS1D and NBINS2D parameters fixed in MonteCarloEngine.h.
    Parameters
    [in]Mod_ithe pointer to the model defined in SomeModel.conf

Definition at line 42 of file MonteCarloEngine.cpp.

42  {
43  TH1D * lhisto = new TH1D("LogLikelihood", "LogLikelihood",
44  NBINS1D, 1., -1.);
45  lhisto->SetDefaultBufferSize(100000);
46  lhisto->GetXaxis()->SetTitle("LogLikelihood");
47  BCH1D * bclhisto = new BCH1D(lhisto);
48  Histo1D["LogLikelihood"] = bclhisto;
49  Mod = Mod_i;
50  int k = 0, kweight = 0;
51  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
52  it < Obs_ALL.end(); it++) {
53  if (!it->isTMCMC()) {
54  k++;
55  if (it->getDistr().compare("noweight") != 0)
56  kweight++;
57  }
58  std::string HistName = it->getName();
59  if (Histo1D.find(HistName) == Histo1D.end()) {
60  TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(),
61  NBINS1D, it->getMin(), it->getMax());
62  histo->GetXaxis()->SetTitle(it->getLabel().c_str());
63  BCH1D * bchisto = new BCH1D(histo);
64  Histo1D[HistName] = bchisto;
65  thMin[it->getName()] = std::numeric_limits<double>::max();
66  thMax[it->getName()] = -std::numeric_limits<double>::max();
67  }
68  }
69  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
70  it < Obs2D_ALL.end(); it++) {
71  if ((it->getDistr()).compare("file") == 0) {
72  if (!it->isTMCMC())
73  throw std::runtime_error("ERROR: cannot handle noMCMC for Observable2D file yet!");
74  } else if (it->getDistr().compare("weight") == 0)
75  throw std::runtime_error("ERROR: do not use Observable2D for analytic 2D weights!");
76  std::string HistName = it->getName();
77  if (Histo2D.find(HistName) == Histo2D.end()) {
78  TH2D * histo2 = new TH2D(HistName.c_str(),
79  (it->getLabel() + " vs " + it->getLabel2()).c_str(),
80  NBINS2D, it->getMin(), it->getMax(),
81  NBINS2D, it->getMin2(), it->getMax2());
82  histo2->GetXaxis()->SetTitle(it->getLabel().c_str());
83  histo2->GetYaxis()->SetTitle(it->getLabel2().c_str());
84  BCH2D * bchisto2 = new BCH2D(histo2);
85  Histo2D[HistName] = bchisto2;
86  }
87  }
88  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
89  it1 != CGO.end(); ++it1) {
90  std::vector<Observable> pino(it1->getObs());
91  for (std::vector<Observable>::iterator it = pino.begin();
92  it != pino.end(); ++it) {
93  //for (int i = 0; i < it1->getObs().size(); i++) {
94  //Observable * it = &(it1->getObs().at(i));
95  if ((it->getDistr()).compare("file") == 0)
96  throw std::runtime_error("Cannot use file in CorrelatedGaussianObservables!");
97  if (!(it->isTMCMC())) {
98  k++;
99  if (it->getDistr().compare("noweight") != 0)
100  throw std::runtime_error("Cannot use weight in CorrelatedGaussianObservables!");
101  }
102  std::string HistName = it->getName();
103  if (Histo1D.find(HistName) == Histo1D.end()) {
104  TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(),
105  NBINS1D, it->getMin(), it->getMax());
106  histo->GetXaxis()->SetTitle(it->getLabel().c_str());
107  BCH1D * bchisto = new BCH1D(histo);
108  Histo1D[HistName] = bchisto;
109  thMin[HistName] = std::numeric_limits<double>::max();
110  thMax[HistName] = -std::numeric_limits<double>::max();
111  }
112  }
113  if (it1->isPrediction()){
114  CorrelationMap[it1->getName()] = new TPrincipal(it1->getObs().size());
115  }
116  }
117  kmax = k;
118  kwmax = kweight;
119 
121 
122 };
std::map< std::string, double > thMax
A map between the name of a theory observable and its maximum value.
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
#define NBINS2D
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
unsigned int kwmax
The number of observables whose weights are used for the MCMC.
std::map< std::string, double > thMin
A map between the name of a theory observable and its minimum value.
void DefineParameters()
A member to classify the prior of the model parameters varied in the Monte Carlo. ...
std::map< std::string, TPrincipal * > CorrelationMap
A map between the name of a theory observable and its maximum value.
std::map< std::string, BCH2D * > Histo2D
A map between pointers to objects of type BCH2D (BAT) and their given names.
unsigned int kmax
The number of observables.
#define NBINS1D
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
Model * Mod
A pointer to an abject of type Model.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
double MonteCarloEngine::LogLikelihood ( const std::vector< double > &  parameters)

This member calculates the loglikelihood for the observables included in the MCMC run.

The model is updated with the new set of parameters through the Model::Update() method. If this is successful this even is counted as a used event. Otherwise it is counted as a discarded event. The log probability is calculated using Weight(const Observable&, const double& ) function for the observables, Weight(const Observable2D&, const double&, const double& ) for the observable2D and Weight(const CorrelatedGaussianObservables& ) for the correlated Gaussian observables. Overloaded from BCEngineMCMC in BAT.

Parameters
[in]parametersthe vector of the parameters that are being varied in the Monte Carlo
Returns
the loglikelihood

Definition at line 248 of file MonteCarloEngine.cpp.

248  {
249  // This methods returns the logarithm of the conditional probability
250  // p(data|parameters). This is where you have to define your model.
251 
252  double logprob = 0.;
253 
254  setDParsFromParameters(parameters, DPars);
255 
256  // if update false set probability equal zero
257  if (!Mod->Update(DPars)) {
258 #ifdef _MCDEBUG
259  std::cout << "event discarded" << std::endl;
260 
261  /* Debug */
262  //for (int k = 0; k < parameters.size(); k++)
263  // std::cout << " " << GetParameter(k)->GetName() << " = "
264  // << DPars[GetParameter(k)->GetName()] << std::endl;
265 #endif
267  return (log(0.));
268  }
269 #ifdef _MCDEBUG
270  //std::cout << "event used in MC" << std::endl;
271 #endif
272 
273  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it != Obs_ALL.end(); it++) {
274  if (it->isTMCMC()) logprob += it->computeWeight();
275  }
276 
277  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it != Obs2D_ALL.end(); it++) {
278  if (it->isTMCMC()) logprob += it->computeWeight();
279  }
280 
281  for (std::vector<CorrelatedGaussianObservables>::iterator it = CGO.begin(); it < CGO.end(); it++) {
282  if(!(it->isPrediction())) logprob += it->computeWeight();
283  }
284  if (isnan(logprob)) {
286  std::cout << "Event discarded since logprob evaluated to NAN.\n";
287  return (log(0.));
288  }
289  NumOfUsedEvents++;
290  return logprob;
291 }
int NumOfDiscardedEvents
The number of events for which the update of the model fails and these events are not used for the MC...
virtual bool Update(const std::map< std::string, double > &DPars)=0
The update method for the model.
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
std::map< std::string, double > DPars
A map between parameter names and their values.
void setDParsFromParameters(const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
int NumOfUsedEvents
The number of events for which the model is successfully updated and hence used for the MCMC run...
complex log(const complex &z)
Model * Mod
A pointer to an abject of type Model.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
void MonteCarloEngine::MCMCIterationInterface ( )

Overloaded from BCEngineMCMC in BAT

The interface is used to update the model parameters using the Model::Update() method. Then the values of the observables are computed and the respective histograms are filled. If the checkTheoryRange flag is true then the minimum and maximum of the theory value is checked and reset to include the current theory values.

Definition at line 293 of file MonteCarloEngine.cpp.

293  {
294 #ifdef _MPI
295  unsigned mychain = 0;
296  int iproc = 0;
297  unsigned npars = GetNParameters();
298  int buffsize = npars + 1;
299  int index_chain[procnum];
300  double *recvbuff = new double[buffsize];
301  std::vector<std::vector<double> > fMCMCxvect;
302  std::vector<double> pars;
303  double **buff;
304 
305 
306  buff = new double*[procnum];
307  int obsbuffsize = Obs_ALL.size() + 2 * Obs2D_ALL.size();
308  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++)
309  obsbuffsize += it1->getObs().size();
310  buff[0] = new double[procnum * obsbuffsize];
311  for (int i = 1; i < procnum; i++) {
312  buff[i] = buff[i - 1] + obsbuffsize;
313  index_chain[i] = -1;
314  }
315 
316  double ** sendbuff = new double *[procnum];
317  sendbuff[0] = new double[procnum * buffsize];
318  for (int il = 1; il < procnum; il++)
319  sendbuff[il] = sendbuff[il - 1] + buffsize;
320 
321  while (mychain < fMCMCNChains) {
322  pars.clear();
323  for (unsigned int k = 0; k < npars; k++)
324  pars.push_back(fMCMCx.at(mychain * npars + k));
325 
326  fMCMCxvect.push_back(pars);
327  index_chain[iproc] = mychain;
328  iproc++;
329  mychain++;
330  if (iproc < procnum && mychain < fMCMCNChains)
331  continue;
332 
333  for (int il = 0; il < iproc; il++) {
334  //The first entry of the array specifies the task to be executed.
335 
336  sendbuff[il][0] = 2.; // 2 = observables calculation
337  for (int im = 1; im < buffsize; im++)
338  sendbuff[il][im] = fMCMCxvect[il][im - 1];
339  }
340  for (int il = iproc; il < procnum; il++) {
341  sendbuff[il][0] = 3.; // 3 = nothing to execute, but return a buffer of observables
342  index_chain[il] = -1;
343  }
344  // double inittime = MPI::Wtime();
345  MPI::COMM_WORLD.Scatter(sendbuff[0], buffsize, MPI::DOUBLE,
346  recvbuff, buffsize, MPI::DOUBLE,
347  0);
348 
349  if (recvbuff[0] == 2.) { // compute observables
350  double sbuff[obsbuffsize];
351  std::map<std::string, double> DPars;
352  pars.assign(recvbuff + 1, recvbuff + buffsize);
353  setDParsFromParameters(pars,DPars);
354  Mod->Update(DPars);
355 
356  int k = 0;
357  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
358  sbuff[k++] = it->computeTheoryValue();
359  }
360  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it < Obs2D_ALL.end(); it++) {
361  sbuff[k++] = it->computeTheoryValue();
362  sbuff[k++] = it->computeTheoryValue2();
363  }
364 
365  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
366  std::vector<Observable> pino(it1->getObs());
367  for (std::vector<Observable>::iterator it = pino.begin(); it != pino.end(); ++it)
368  sbuff[k++] = it->computeTheoryValue();
369  }
370  MPI::COMM_WORLD.Gather(sbuff, obsbuffsize, MPI::DOUBLE, buff[0], obsbuffsize, MPI::DOUBLE, 0);
371  } else if (recvbuff[0] == 3.) { // do not compute observables, but gather the buffer
372  double sbuff[obsbuffsize];
373  MPI::COMM_WORLD.Gather(sbuff, obsbuffsize, MPI::DOUBLE, buff[0], obsbuffsize, MPI::DOUBLE, 0);
374  }
375 
376  for (int il = 0; il < procnum; il++) {
377  if (index_chain[il] >= 0) {
378  int k = 0;
379  // fill the histograms for observables
380  int ko = 0, kweight = 0;
381  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
382  it < Obs_ALL.end(); it++) {
383  double th = buff[il][k++];
384  /* set the min and max of theory values */
385  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
386  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
387  Histo1D[it->getName()]->GetHistogram()->Fill(th);
388  if (!it->isTMCMC()) {
389  obval[index_chain[il] * kmax + ko++] = th;
390  if (it->getDistr().compare("noweight") != 0) {
391  obweight[index_chain[il] * kwmax + kweight++] = it->computeWeight(th);
392  }
393  }
394  }
395 
396  // fill the 2D histograms for observables
397  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
398  it < Obs2D_ALL.end(); it++) {
399  double th1 = buff[il][k++];
400  double th2 = buff[il][k++];
401  Histo2D[it->getName()]->GetHistogram()->Fill(th1, th2);
402  }
403 
404  // fill the histograms for correlated observables
405  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
406  it1 < CGO.end(); it1++) {
407  std::vector<Observable> pino(it1->getObs());
408  Double_t * COdata = new Double_t[pino.size()];
409  int nObs = 0;
410  for (std::vector<Observable>::iterator it = pino.begin();
411  it != pino.end(); ++it) {
412  double th = buff[il][k++];
413  /* set the min and max of theory values */
414  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
415  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
416  Histo1D[it->getName()]->GetHistogram()->Fill(th);
417  if (it1->isPrediction()) COdata[nObs++] = th;
418  }
419  if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
420  delete [] COdata;
421  }
422  }
423  }
424  iproc = 0;
425  fMCMCxvect.clear();
426  }
427  delete sendbuff[0];
428  delete [] sendbuff;
429  delete [] recvbuff;
430  delete buff[0];
431  delete [] buff;
432 #else
433  for (unsigned int i = 0; i < fMCMCNChains; ++i) {
434  std::vector<double>::const_iterator first = fMCMCx.begin() + i * GetNParameters();
435  std::vector<double>::const_iterator last = first + GetNParameters();
436  std::vector<double> currvec(first, last);
437  setDParsFromParameters(currvec,DPars);
438 
439  Mod->Update(DPars);
440 
441  // fill the histograms for observables
442  int k = 0, kweight = 0;
443  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
444  it < Obs_ALL.end(); it++) {
445  double th = it->computeTheoryValue();
446  /* set the min and max of theory values */
447  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
448  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
449  Histo1D[it->getName()]->GetHistogram()->Fill(th);
450  if (!it->isTMCMC()) {
451  obval[i * kmax + k] = th;
452  k++;
453  if (it->getDistr().compare("noweight") != 0) {
454  obweight[i * kwmax + kweight] = it->computeWeight();
455  kweight++;
456  }
457  }
458  }
459 
460  // fill the 2D histograms for observables
461  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
462  it < Obs2D_ALL.end(); it++) {
463  double th1 = it->computeTheoryValue();
464  double th2 = it->computeTheoryValue2();
465  Histo2D[it->getName()]->GetHistogram()->Fill(th1, th2);
466  }
467 
468  // fill the histograms for correlated observables
469  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
470  it1 < CGO.end(); it1++) {
471  std::vector<Observable> pino(it1->getObs());
472  Double_t * COdata = new Double_t[pino.size()];
473  int nObs = 0;
474  for (std::vector<Observable>::iterator it = pino.begin();
475  it != pino.end(); ++it) {
476  double th = it->computeTheoryValue();
477  /* set the min and max of theory values */
478  if (th < thMin[it->getName()]) thMin[it->getName()] = th;
479  if (th > thMax[it->getName()]) thMax[it->getName()] = th;
480  Histo1D[it->getName()]->GetHistogram()->Fill(th);
481  if (it1->isPrediction()) COdata[nObs++] = th;
482  }
483  if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
484  delete [] COdata;
485  }
486  }
487 #endif
488  for (unsigned int i = 0; i < fMCMCNChains; i++)
489  {
490  Histo1D["LogLikelihood"]->GetHistogram()->Fill(MCMCGetLogProbx(i)-LogAPrioriProbability(MCMCGetx(i)));
491  }
492 }
virtual bool Update(const std::map< std::string, double > &DPars)=0
The update method for the model.
std::map< std::string, double > thMax
A map between the name of a theory observable and its maximum value.
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
unsigned int kwmax
The number of observables whose weights are used for the MCMC.
std::map< std::string, double > thMin
A map between the name of a theory observable and its minimum value.
double * obval
A pointer to the vector of observable values.
std::map< std::string, double > DPars
A map between parameter names and their values.
void setDParsFromParameters(const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
std::map< std::string, TPrincipal * > CorrelationMap
A map between the name of a theory observable and its maximum value.
std::map< std::string, BCH2D * > Histo2D
A map between pointers to objects of type BCH2D (BAT) and their given names.
double * obweight
A pointer to the vector of observable weights.
unsigned int kmax
The number of observables.
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
Model * Mod
A pointer to an abject of type Model.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
void MonteCarloEngine::PrintCorrelationMatrix ( const std::string  filename)

This member generates the correlation matrix using BCH2D from the BAT libraries.

Parameters
[in]filenamethe name of the file where the correlation matrix is printed

Definition at line 611 of file MonteCarloEngine.cpp.

611  {
612  std::ofstream out;
613  out.open(filename.c_str(), std::ios::out);
614 
615  int npar = GetNParameters();
616 
617  for (int i = 0; i < npar; ++i)
618  out << " & " << GetParameter(i)->GetName();
619  out << " \\\\" << std::endl;
620 
621  for (int i = 0; i < npar; ++i) {
622  out << GetParameter(i)->GetName() << " & $";
623  for (int j = 0; j < npar; ++j) {
624  if (i != j) {
625  BCH2D* bch2d_temp = GetMarginalized(GetParameter(i), GetParameter(j));
626  if (bch2d_temp != NULL)
627  out << bch2d_temp->GetHistogram()->GetCorrelationFactor();
628  else
629  out << 0.;
630  } else
631  out << 1.;
632  if (j == npar - 1) out << "$ \\\\" << std::endl;
633  else out << "$ & $";
634  }
635  }
636 
637  out.close();
638 }
void MonteCarloEngine::PrintHistogram ( BCModelOutput &  out,
Observable it,
const std::string  OutputDir 
)

Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables.

Parameters
[in]outa reference to an object of type BCModelOutput as defined in the BAT libraries
[in]ita reference to an object of type Observable
[in]OutputDirthe name of the output directory

Definition at line 520 of file MonteCarloEngine.cpp.

521  {
522  std::string HistName = it.getName();
523  double min = thMin[it.getName()];
524  double max = thMax[it.getName()];
525 
526  if (Histo1D[HistName]->GetHistogram()->Integral() > 0.0) {
527  std::string fname = OutputDir + "/" + HistName + ".pdf";
528  // BCH1D* pippo = Histo1D[HistName];
529  // double x = pippo->GetMean();
530  // pippo->Print("Dmd1.pdf");
531  Histo1D[HistName]->SetGlobalMode(it.computeTheoryValue());
532  Histo1D[HistName]->Print(fname.c_str());
533  std::cout << fname << " has been created." << std::endl;
534  out.Write(Histo1D[HistName]->GetHistogram());
535  CheckHistogram(*Histo1D[HistName]->GetHistogram(), it.getName());
536  } else
537  HistoLog << "WARNING: The histogram of "
538  << it.getName() << " is empty!" << std::endl;
539 
540  HistoLog.precision(10);
541  HistoLog << " [min, max]=[" << min << ", " << max << "]" << std::endl;
542  HistoLog.precision(6);
543 }
std::map< std::string, double > thMax
A map between the name of a theory observable and its maximum value.
std::map< std::string, double > thMin
A map between the name of a theory observable and its minimum value.
void CheckHistogram(const TH1D &hist, const std::string name)
This member checks if there is overflow of the 1D histogram.
std::string getName() const
A get method to access the name of the observable.
Definition: Observable.h:304
double computeTheoryValue()
A method to access the computed theory value of the observable.
Definition: Observable.cpp:115
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
void MonteCarloEngine::PrintHistogram ( BCModelOutput &  out,
const std::string  OutputDir 
)

Member used for printing histograms for observables, observable2D, correlated Gaussian observables and model parameters vs. observables.

Parameters
[in]outa reference to an object of type BCModelOutput as defined in the BAT libraries
[in]OutputDirthe name of the output directory

Definition at line 545 of file MonteCarloEngine.cpp.

545  {
546  std::vector<double> mode(GetBestFitParameters());
547  if (mode.size() == 0) {
548  if (rank == 0) throw std::runtime_error("\n ERROR: Global Mode could not be determined possibly because of infinite loglikelihood. Observables histogram cannot be generated.\n");
549  else sleep(2);
550  }
552 
553  Mod->Update(DPars);
554 
555  // print the histograms to pdf files
556  for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end();
557  it++) {
558  PrintHistogram(out, *it, OutputDir);
559  }
560  for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
561  it1 < CGO.end(); it1++) {
562  std::vector<Observable> pino(it1->getObs());
563  for (std::vector<Observable>::iterator it = pino.begin();
564  it != pino.end(); ++it)
565  PrintHistogram(out, *it, OutputDir);
566  }
567  for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
568  it < Obs2D_ALL.end(); it++) {
569  std::string HistName = it->getName();
570  if (Histo2D[HistName]->GetHistogram()->Integral() > 0.0) {
571  std::string fname = OutputDir + "/" + HistName + ".pdf";
572  double th[2];
573  th[0] = it->computeTheoryValue();
574  th[1] = it->computeTheoryValue2();
575  Histo2D[HistName]->SetGlobalMode(th);
576  Histo2D[HistName]->Print(fname.c_str());
577  std::cout << fname << " has been created." << std::endl;
578  out.Write(Histo2D[HistName]->GetHistogram());
579  CheckHistogram(*Histo2D[HistName]->GetHistogram(), HistName);
580  } else
581  HistoLog << "WARNING: The histogram of "
582  << HistName << " is empty!" << std::endl;
583  }
584  std::string fname = OutputDir + "/LogLikelihood.pdf";
585  Histo1D["LogLikelihood"]->Print(fname.c_str());
586  std::cout << fname << " has been created." << std::endl;
587  out.Write(Histo1D["LogLikelihood"]->GetHistogram());
588 }
virtual bool Update(const std::map< std::string, double > &DPars)=0
The update method for the model.
std::vector< Observable2D > & Obs2D_ALL
A vector of all pairs of observable for Observable2D.
boost::ptr_vector< Observable > & Obs_ALL
A vector of all observables.
void CheckHistogram(const TH1D &hist, const std::string name)
This member checks if there is overflow of the 1D histogram.
std::map< std::string, double > DPars
A map between parameter names and their values.
void PrintHistogram(BCModelOutput &out, Observable &it, const std::string OutputDir)
Overloaded from PrintHistogram(BCModelOutput&, const std::string) to print histogram for observables...
void setDParsFromParameters(const std::vector< double > &parameters, std::map< std::string, double > &DPars_i)
int rank
Rank of the process for a MPI run. Value is 0 for a serial run.
std::map< std::string, BCH2D * > Histo2D
A map between pointers to objects of type BCH2D (BAT) and their given names.
std::map< std::string, BCH1D * > Histo1D
A map between pointers to objects of type BCH1D (BAT) and their given names.
std::ostringstream HistoLog
A stream to store the output messages from printing and checking histograms.
Model * Mod
A pointer to an abject of type Model.
std::vector< CorrelatedGaussianObservables > & CGO
A vector of correlated Gaussian observables.
double MonteCarloEngine::SecondDerivative ( const BCParameter *  par1,
const BCParameter *  par2,
std::vector< double >  point 
)

Definition at line 849 of file MonteCarloEngine.cpp.

849  {
850 
851  if (point.size() != GetNParameters()) {
852  throw std::runtime_error("MCMCENgine::SecondDerivative : Invalid number of entries in the vector.");
853  }
854 
855  // define steps
856  const double dy1 = par2->GetRangeWidth() / NSTEPS;
857  const double dy2 = dy1 * 2.;
858  const double dy3 = dy1 * 3.;
859 
860  // define points at which to evaluate
861  std::vector<double> y1p = point;
862  std::vector<double> y1m = point;
863  std::vector<double> y2p = point;
864  std::vector<double> y2m = point;
865  std::vector<double> y3p = point;
866  std::vector<double> y3m = point;
867 
868  unsigned idy = GetParameters().Index(par2->GetName());
869 
870  y1p[idy] += dy1;
871  y1m[idy] -= dy1;
872  y2p[idy] += dy2;
873  y2m[idy] -= dy2;
874  y3p[idy] += dy3;
875  y3m[idy] -= dy3;
876 
877  const double m1 = (FirstDerivative(par1, y1p) - FirstDerivative(par1, y1m)) / 2. / dy1;
878  const double m2 = (FirstDerivative(par1, y2p) - FirstDerivative(par1, y2m)) / 4. / dy1;
879  const double m3 = (FirstDerivative(par1, y3p) - FirstDerivative(par1, y3m)) / 6. / dy1;
880 
881  return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
882 }
double FirstDerivative(const BCParameter *par, std::vector< double > point)
#define NSTEPS
void MonteCarloEngine::setDParsFromParameters ( const std::vector< double > &  parameters,
std::map< std::string, double > &  DPars_i 
)

Definition at line 199 of file MonteCarloEngine.cpp.

201 {
202  std::map<std::string, std::vector<double> > cgpmap;
203 
204  unsigned int k = 0;
205  for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it < ModPars.end(); it++){
206  if(it->IsFixed())
207  continue;
208  if(it->getname().compare(GetParameter(k)->GetName()) != 0)
209  {
210  std::stringstream out;
211  out << it->getname();
212  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + "is sitting at the wrong position in the BAT parameters vector");
213  }
214  if (it->IsCorrelated()) {
215  std::string index = it->getname().substr(it->getCgp_name().size());
216  unsigned long int lindex = strtol(index.c_str(),NULL,10);
217  if (lindex - 1 == cgpmap[it->getCgp_name()].size())
218  cgpmap[it->getCgp_name()].push_back(parameters[k]);
219  else {
220  std::stringstream out;
221  out << it->getname() << " " << lindex;
222  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + "seems to be a CorrelatedGaussianParameters object but the corresponding parameters are missing or not in the right order");
223  }
224 
225  } else
226  DPars_i[it->getname()] = parameters[k];
227  k++;
228  }
229 
230  for (unsigned int j = 0; j < CGP.size(); j++) {
231  std::vector<double> current = cgpmap.at(CGP[j].getName());
232  if (current.size() != CGP[j].getPars().size()) {
233  std::stringstream out;
234  out << CGP[j].getName();
235  throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + " appears to be represented in cgpmap with a wrong size");
236  }
237 
238  std::vector<double> porig = CGP[j].getOrigParsValue(current);
239 
240  for(unsigned int l = 0; l < porig.size(); l++) {
241  DPars_i[CGP[j].getPar(l).getname()] = porig[l];
242  }
243  }
244 }
const std::vector< ModelParameter > & ModPars
A vector of model parameters.
const std::vector< CorrelatedGaussianParameters > & CGP
A vector of correlated Gaussian parameters.
void MonteCarloEngine::setHistogramOverFlow ( bool  overflow)
inline

Definition at line 235 of file MonteCarloEngine.h.

236  {
237  if (overflow) TH1::StatOverflows(kTRUE);
238  else TH1::StatOverflows(kFALSE);
239  }
void MonteCarloEngine::setNChains ( unsigned int  i)

A set method to fix the number of chains.

The number of chains are set using the MCMCSetNChains() from BCEngineMCMC class in BAT. This also creates pointers to the vector of observable values and weights for all the chains.

Parameters
[in]ithe number of chains

Definition at line 124 of file MonteCarloEngine.cpp.

124  {
125  MCMCSetNChains(i);
126  obval = new double[fMCMCNChains * kmax];
127  obweight = new double[fMCMCNChains * kwmax];
128 }
unsigned int kwmax
The number of observables whose weights are used for the MCMC.
double * obval
A pointer to the vector of observable values.
double * obweight
A pointer to the vector of observable weights.
unsigned int kmax
The number of observables.
std::string MonteCarloEngine::writePreRunData ( )

A method to write in a text file the best fit parameters and the prerun scale factors.

Returns
a string containing the data

Definition at line 822 of file MonteCarloEngine.cpp.

823 {
824  std::vector<double> mode(GetBestFitParameters());
825  std::vector<double> scales(MCMCGetTrialFunctionScaleFactor(0));
826  std::ostringstream StatsLog;
827  for (unsigned int i = 0; i < mode.size(); i++)
828  StatsLog << GetParameter(i)->GetName() << " " << mode.at(i) << " " << scales.at(i) << std::endl;
829  return StatsLog.str().c_str();
830 }

Member Data Documentation

std::vector<CorrelatedGaussianObservables>& MonteCarloEngine::CGO
private

A vector of correlated Gaussian observables.

Definition at line 256 of file MonteCarloEngine.h.

const std::vector<CorrelatedGaussianParameters>& MonteCarloEngine::CGP
private

A vector of correlated Gaussian parameters.

Definition at line 253 of file MonteCarloEngine.h.

std::map<std::string, TPrincipal *> MonteCarloEngine::CorrelationMap
private

A map between the name of a theory observable and its maximum value.

Definition at line 263 of file MonteCarloEngine.h.

std::map<std::string, double> MonteCarloEngine::DPars
private

A map between parameter names and their values.

Definition at line 258 of file MonteCarloEngine.h.

std::map<std::string, BCH1D * > MonteCarloEngine::Histo1D
private

A map between pointers to objects of type BCH1D (BAT) and their given names.

Definition at line 259 of file MonteCarloEngine.h.

std::map<std::string, BCH2D * > MonteCarloEngine::Histo2D
private

A map between pointers to objects of type BCH2D (BAT) and their given names.

Definition at line 260 of file MonteCarloEngine.h.

std::ostringstream MonteCarloEngine::HistoLog
private

A stream to store the output messages from printing and checking histograms.

Definition at line 268 of file MonteCarloEngine.h.

unsigned int MonteCarloEngine::kmax
private

The number of observables.

Definition at line 267 of file MonteCarloEngine.h.

unsigned int MonteCarloEngine::kwmax
private

The number of observables whose weights are used for the MCMC.

Definition at line 266 of file MonteCarloEngine.h.

Model* MonteCarloEngine::Mod
private

A pointer to an abject of type Model.

Definition at line 257 of file MonteCarloEngine.h.

const std::vector<ModelParameter>& MonteCarloEngine::ModPars
private

A vector of model parameters.

Definition at line 252 of file MonteCarloEngine.h.

int MonteCarloEngine::NumOfDiscardedEvents
private

The number of events for which the update of the model fails and these events are not used for the MCMC run.

Definition at line 270 of file MonteCarloEngine.h.

int MonteCarloEngine::NumOfUsedEvents
private

The number of events for which the model is successfully updated and hence used for the MCMC run.

Definition at line 269 of file MonteCarloEngine.h.

std::vector<Observable2D>& MonteCarloEngine::Obs2D_ALL
private

A vector of all pairs of observable for Observable2D.

Definition at line 255 of file MonteCarloEngine.h.

boost::ptr_vector<Observable>& MonteCarloEngine::Obs_ALL
private

A vector of all observables.

Definition at line 254 of file MonteCarloEngine.h.

double* MonteCarloEngine::obval
private

A pointer to the vector of observable values.

Definition at line 264 of file MonteCarloEngine.h.

double* MonteCarloEngine::obweight
private

A pointer to the vector of observable weights.

Definition at line 265 of file MonteCarloEngine.h.

int MonteCarloEngine::rank
private

Rank of the process for a MPI run. Value is 0 for a serial run.

Definition at line 271 of file MonteCarloEngine.h.

std::map<std::string, double> MonteCarloEngine::thMax
private

A map between the name of a theory observable and its maximum value.

Definition at line 262 of file MonteCarloEngine.h.

std::map<std::string, double> MonteCarloEngine::thMin
private

A map between the name of a theory observable and its minimum value.

Definition at line 261 of file MonteCarloEngine.h.


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