MonteCarloEngine Class Reference

An engine class for Monte Carlo. More...

#include <MonteCarloEngine.h>

Inheritance diagram for MonteCarloEngine:

Collaboration diagram for MonteCarloEngine:

Detailed Description

An engine class for Monte Carlo.

Author

HEPfit Collaboration

Copyright

GNU General Public License

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...
Model *	Mod
	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_i	the vector of model parameters as defined in SomeModel.conf
[in]	Obs_i	the vector of observables as defined in SomeModel.conf
[in]	Obs2D_i	the vector of observables2D as defined in SomeModel.conf
[in]	CGO_i	the vector of correlated Gaussian observables as defined in SomeModel.conf

Definition at line 23 of file MonteCarloEngine.cpp.

: BCModel(""), ModPars(ModPars_i), CGP(CGP_i), Obs_ALL(Obs_i), Obs2D_ALL(Obs2D_i),
  CGO(CGO_i), NumOfUsedEvents(0), NumOfDiscardedEvents(0) {
    obval = NULL;
    obweight = NULL;
    Mod = NULL;
    TH1::StatOverflows(kTRUE);
#ifdef _MPI
    rank = MPI::COMM_WORLD.Get_rank();
#else
    rank = 0;
#endif
};

MonteCarloEngine::~MonteCarloEngine

(

)

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

Definition at line 132 of file MonteCarloEngine.cpp.

{
    delete [] obval;
    delete [] obweight;
    /* The following code has been commented out pending further review.
       It is causing crashes at the termination of the code if the histograms
       are accessed from the main program.*/
    //    for (std::map<std::string, BCH1D *>::iterator it = Histo1D.begin();
    //            it != Histo1D.end(); it++)
    //        delete it->second;
    //    for (std::map<std::string, BCH2D *>::iterator it = Histo2D.begin();
    //            it != Histo2D.end(); it++)
    //        delete it->second;
    for (std::map<std::string, TPrincipal *>::iterator it = CorrelationMap.begin();
            it != CorrelationMap.end(); it++)
        delete it->second;
};

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.

                                 {
    fMCMCFlagWritePreRunToFile = false;
    int k = 0, kweight = 0;
    for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
            it < Obs_ALL.end(); it++) {
        if (!it->isTMCMC()) {
            for (unsigned int i = 0; i < fMCMCNChains; ++i)
                fMCMCTrees[i]->Branch(it->getName().c_str(), &obval[i * kmax + k],
                    (it->getName() + "/D").c_str());
            k++;
            if (it->getDistr().compare("noweight") != 0) {
                for (unsigned int i = 0; i < fMCMCNChains; ++i)
                    fMCMCTrees[i]->Branch((it->getName() + "_weight").c_str(),
                        &obweight[i * kwmax + kweight],
                        (it->getName() + "_weight/D").c_str());
                kweight++;
            }
        }
    }
}

void MonteCarloEngine::CheckHistogram	(	const TH1D &	hist,
		const std::string	name
	)

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

Parameters

[in]	hist	a reference to an object of type TH1D as defined in the ROOT libraries
[in]	name	the name for the histogram

Definition at line 494 of file MonteCarloEngine.cpp.

                                                                            {
    // output the portions of underflow and overflow bins
    double UnderFlowContent = hist.GetBinContent(0);
    double OverFlowContent = hist.GetBinContent(NBINS1D + 1);
    double Integral = hist.Integral();
    double TotalContent = 0.0;
    for (int n = 0; n <= NBINS1D + 1; n++)
        TotalContent += hist.GetBinContent(n);
    HistoLog << name << ": "
            << Integral / TotalContent * 100. << "% within the range, "
            << UnderFlowContent / TotalContent * 100. << "% underflow, "
            << OverFlowContent / TotalContent * 100. << "% overflow"
            << std::endl;
}

void MonteCarloEngine::CheckHistogram	(	const TH2D &	hist,
		const std::string	name
	)

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

Parameters

[in]	hist	a reference to an object of type TH2D as defined in the ROOT libraries
[in]	name	the name for the histogram

Definition at line 509 of file MonteCarloEngine.cpp.

                                                                            {
    double Integral = hist.Integral();
    double TotalContent = 0.0;
    for (int m = 0; m <= NBINS2D + 1; m++)
        for (int n = 0; n <= NBINS2D + 1; n++)
            TotalContent += hist.GetBinContent(m, n);
    HistoLog << name << ": "
            << Integral / TotalContent * 100. << "% within the ranges"
            << std::endl;
}

double MonteCarloEngine::computeNormalization

(

)

Definition at line 832 of file MonteCarloEngine.cpp.

                                              {

    unsigned int Npars = GetNParameters();
    std::vector<double> mode(GetBestFitParameters());
    gslpp::matrix<double> Hessian(Npars, Npars, 0.);

    for (unsigned int i = 0; i < Npars; i++)
        for (unsigned int j = 0; j < Npars; j++) {
            // calculate Hessian matrix element
            Hessian.assign(i, j, -SecondDerivative(GetParameter(i), GetParameter(j), GetBestFitParameters()));
        }

    double det_Hessian = Hessian.determinant();

    return exp(Npars / 2. * log(2. * M_PI) + 0.5 * log(1. / det_Hessian) + LogLikelihood(mode) + LogAPrioriProbability(mode));
}

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.

                                              {
    std::vector<double> mode(GetBestFitParameters());
    if (mode.size() == 0) {
        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");
        else sleep (2);
    }
    std::ostringstream StatsLog;
    int i = 0;
    StatsLog << "Statistics file for Observables, Binned Observables and Corellated Gaussian Observables.\n" << std::endl;
    if (Obs_ALL.size() > 0) StatsLog << "Observables:\n" << std::endl;
    for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {

        if (it->getObsType().compare("BinnedObservable") == 0) {
            StatsLog << "  (" << ++i << ") Binned Observable \"";
            StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << ", " << it->getTho()->getBinMax() << "]" << "\":";
        } else if (it->getObsType().compare("FunctionObservable") == 0) {
            StatsLog << "  (" << ++i << ") Function Observable \"";
            StatsLog << it->getName() << "[" << it->getTho()->getBinMin() << "]" << "\":";
        } else if (it->getObsType().compare("HiggsObservable") == 0) {
            StatsLog << "  (" << ++i << ") Higgs Observable \"";
            StatsLog << it->getName() << "\":";
        } else {
            StatsLog << "  (" << ++i << ") " << it->getObsType() << " \"";
            StatsLog << it->getName() << "\":";
        }

        StatsLog << std::endl;

        BCH1D * bch1d = Histo1D[it->getName()];
        if (bch1d->GetHistogram()->Integral() > 0.0) {
            StatsLog << "      Mean +- sqrt(V):                " << std::setprecision(4)
                    << bch1d->GetMean() << " +- " << std::setprecision(4)
                    << bch1d->GetRMS() << std::endl

                    << "      (Marginalized) mode:            " << bch1d->GetMode() << std::endl;

            std::vector<double> v1;
            std::vector<double> v2;
            std::vector<double> v3;
            v1 = bch1d->GetSmallestIntervals(0.6827);
            v2 = bch1d->GetSmallestIntervals(0.9545);
            v3 = bch1d->GetSmallestIntervals(0.9973);
            StatsLog << "      Smallest interval(s) containing at least 68.27% and local mode(s):"
                    << std::endl;
            for (unsigned j = 0; j < v1.size(); j += 5)
                StatsLog << "       (" << v1[j] << ", " << v1[j + 1]
                    << ") (local mode at " << v1[j + 3] << " with rel. height "
                    << v1[j + 2] << "; rel. area " << v1[j + 4] << ")"
                    << std::endl;
            StatsLog << std::endl;

            StatsLog << "      Smallest interval(s) containing at least 95.45% and local mode(s):"
                    << std::endl;
            for (unsigned j = 0; j < v2.size(); j += 5)
                StatsLog << "       (" << v2[j] << ", " << v2[j + 1]
                    << ") (local mode at " << v2[j + 3] << " with rel. height "
                    << v2[j + 2] << "; rel. area " << v2[j + 4] << ")"
                    << std::endl;
            StatsLog << std::endl;

            StatsLog << "      Smallest interval(s) containing at least 99.73% and local mode(s):"
                    << std::endl;
            for (unsigned j = 0; j < v3.size(); j += 5)
                StatsLog << "       (" << v3[j] << ", " << v3[j + 1]
                    << ") (local mode at " << v3[j + 3] << " with rel. height "
                    << v3[j + 2] << "; rel. area " << v3[j + 4] << ")"
                    << std::endl;
            StatsLog << std::endl;
        } else {
            StatsLog << "\nWARNING: The histogram of " << it->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
        }
    }

    if (CGO.size() > 0) StatsLog << "\nCorellated (Gaussian) Observables:\n" << std::endl;
    for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
            it1 < CGO.end(); it1++) {
        StatsLog << "\n" << it1->getName() << ":\n" << std::endl;
        i = 0;
        std::vector<Observable> CGObs(it1->getObs());
        for (std::vector<Observable>::iterator it2 = CGObs.begin(); it2 < CGObs.end(); it2++) {

            if (it2->getObsType().compare("BinnedObservable") == 0) {
                StatsLog << "  (" << ++i << ") Binned Observable \"";
                StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << ", " << it2->getTho()->getBinMax() << "]" << "\":";
            } else if (it2->getObsType().compare("FunctionObservable") == 0) {
                StatsLog << "  (" << ++i << ") Function Observable \"";
                StatsLog << it2->getName() << "[" << it2->getTho()->getBinMin() << "]" << "\":";
            } else if (it2->getObsType().compare("HiggsObservable") == 0) {
                StatsLog << "  (" << ++i << ") Higgs Observable \"";
                StatsLog << it2->getName() << "\":";
            } else {
                StatsLog << "  (" << ++i << ") " << it2->getObsType() << " \"";
                StatsLog << it2->getName() << "\":";
            }

            StatsLog << std::endl;
            BCH1D * bch1d = Histo1D[it2->getName()];
            if (bch1d->GetHistogram()->Integral() > 0.0) {
                StatsLog << "      Mean +- sqrt(V):                " << std::setprecision(4)
                        << bch1d->GetMean() << " +- " << std::setprecision(4)
                        << bch1d->GetRMS() << std::endl

                        << "      (Marginalized) mode:            " << bch1d->GetMode() << std::endl;

                std::vector<double> v1;
                std::vector<double> v2;
                std::vector<double> v3;
                v1 = bch1d->GetSmallestIntervals(0.6827);
                v2 = bch1d->GetSmallestIntervals(0.9545);
                v3 = bch1d->GetSmallestIntervals(0.9973);
                StatsLog << "      Smallest interval(s) containing at least 68.27% and local mode(s):"
                        << std::endl;
                for (unsigned j = 0; j < v1.size(); j += 5)
                    StatsLog << "       (" << v1[j] << ", " << v1[j + 1]
                        << ") (local mode at " << v1[j + 3] << " with rel. height "
                        << v1[j + 2] << "; rel. area " << v1[j + 4] << ")"
                        << std::endl;
                StatsLog << std::endl;

                StatsLog << "      Smallest interval(s) containing at least 95.45% and local mode(s):"
                        << std::endl;
                for (unsigned j = 0; j < v2.size(); j += 5)
                    StatsLog << "       (" << v2[j] << ", " << v2[j + 1]
                        << ") (local mode at " << v2[j + 3] << " with rel. height "
                        << v2[j + 2] << "; rel. area " << v2[j + 4] << ")"
                        << std::endl;
                StatsLog << std::endl;

                StatsLog << "      Smallest interval(s) containing at least 99.73% and local mode(s):"
                        << std::endl;
                for (unsigned j = 0; j < v3.size(); j += 5)
                    StatsLog << "       (" << v3[j] << ", " << v3[j + 1]
                        << ") (local mode at " << v3[j + 3] << " with rel. height "
                        << v3[j + 2] << "; rel. area " << v3[j + 4] << ")"
                        << std::endl;
                StatsLog << std::endl;
            } else {
                StatsLog << "\nWARNING: The histogram of " << it2->getName() << " is empty! Statistics cannot be generated\n" << std::endl;
            }
        }
        if (it1->isPrediction()) {
            int size = it1->getObs().size();
            CorrelationMap[it1->getName()]->MakePrincipals();
            //CorrelationMap[it1->getName()]->Print();
            TMatrixD * corr = const_cast<TMatrixD*>(CorrelationMap[it1->getName()]->GetCovarianceMatrix());
            *corr *= (double)size;
            StatsLog << "\nThe correlation matrix for " << it1->getName() << " is given by the " << size << "x"<< size << " matrix:\n" << std::endl;

            for (int i = 0; i < size + 1; i++) {
                if (i == 0) StatsLog << std::setw(4) << "" << " | ";
                else StatsLog << std::setw(5) << i << std::setw(5) << "    |";
            }
            StatsLog << std::endl;
            for (int i = 0; i < size + 1; i++) {
                if (i == 0) StatsLog << std::setw(8) << "--------";
                else StatsLog << std::setw(10) << "----------";
            }
            StatsLog << std::endl;
            for (int i = 0; i < size; i++) {
                for (int j = 0; j < size + 1; j++) {
                    if (j == 0) StatsLog << std::setw(4) << i+1 << " |";
                    else StatsLog << std::setprecision(5) << std::setw(10) << (*corr)(i, j - 1);
                }
            StatsLog << std::endl;
            }
        }
    }

    std::vector<double> parsa;
    for (unsigned int i = 0; i < GetNParameters(); i++)
        parsa.push_back(GetMarginalized(fParameters[i])->GetMean());
    double llik = LogLikelihood(parsa);
    StatsLog << "LogLikelihood on the mean values of parameters: " << llik << std::endl;
    double llika = Histo1D["LogLikelihood"]->GetMean();
    StatsLog << "LogLikelihood mean value: " << llika << std::endl;
    double dbar = -2.*llika; //Wikipedia notation... 
    double dothetabar = -2.*llik; //Wikipedia notation...
    StatsLog << "IC value (don't ask me what it means...): " << 3.*dbar - 2.*dothetabar << std::endl; 
    StatsLog << "DIC value (same as above...): " << 2.*dbar - dothetabar << std::endl; 
    return StatsLog.str().c_str();
}

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.

                                        {
    // Add parameters to your model here.
    // You can then use them in the methods below by calling the
    // parameters.at(i) or parameters[i], where i is the index
    // of the parameter. The indices increase from 0 according to the
    // order of adding the parameters.
    if (rank == 0) std::cout << "\nParameters varied in Monte Carlo:" << std::endl;
    int k = 0;
    for (std::vector<ModelParameter>::const_iterator it = ModPars.begin();
            it < ModPars.end(); it++) {
        if (it->geterrf() == 0. && it->geterrg() == 0.)
            continue;
        AddParameter(it->getname().c_str(), it->getmin(), it->getmax());
        if (rank == 0) std::cout << k << ": " << it->getname() << ", ";
        if (it->IsCorrelated()) {
            for (unsigned int i = 0; i < CGP.size(); i++) {
                if (CGP[i].getName().compare(it->getCgp_name()) == 0) {
                    std::string index = it->getname().substr(CGP[i].getName().size());
                    long int lindex = strtol(index.c_str(),NULL,10);
                    if (lindex > 0)
                        DPars[CGP[i].getPar(lindex - 1).getname()] = 0.;
                    else {
                        std::stringstream out;
                        out << it->getname();
                        throw std::runtime_error("MonteCarloEngine::DefineParameters(): " + out.str() + "seems to be part of a CorrelatedGaussianParameters object, but I couldn't find the corresponding object");
                    }
                }
            }
        } else
            DPars[it->getname()] = 0.;
        if (it->geterrf() == 0.) SetPriorGauss(k, it->getave(), it->geterrg());
        else if (it->geterrg() == 0.) SetPriorConstant(k);
        else {
            TF1 * combined = new TF1(it->getname().c_str(),
                    "(TMath::Erf((x-[0]+[2])/sqrt(2.)/[1])-TMath::Erf((x-[0]-[2])/sqrt(2.)/[1]))/4./[2]",
                    it->getmin(), it->getmax());
            combined->SetParameter(0, it->getave());
            combined->SetParameter(1, it->geterrg());
            combined->SetParameter(2, it->geterrf());
            SetPrior(k, combined);
            delete combined;
        }
        k++;
    }
}

double MonteCarloEngine::FirstDerivative	(	const BCParameter *	par,
		std::vector< double >	point
	)

Definition at line 884 of file MonteCarloEngine.cpp.

                                                                                         {

    if (point.size() != GetNParameters()) {
        throw std::runtime_error("MCMCENgine::FirstDerivative : Invalid number of entries in the vector.");
    }

    // define steps
    const double dx1 = par->GetRangeWidth() / NSTEPS;
    const double dx2 = dx1 * 2.;
    const double dx3 = dx1 * 3.;

    // define points at which to evaluate
    std::vector<double> x1p = point;
    std::vector<double> x1m = point;
    std::vector<double> x2p = point;
    std::vector<double> x2m = point;
    std::vector<double> x3p = point;
    std::vector<double> x3m = point;

    unsigned idx = GetParameters().Index(par->GetName());

    x1p[idx] += dx1;
    x1m[idx] -= dx1;
    x2p[idx] += dx2;
    x2m[idx] -= dx2;
    x3p[idx] += dx3;
    x3m[idx] -= dx3;

    const double m1 = (Function_h(x1p) - Function_h(x1m)) / 2. / dx1;
    const double m2 = (Function_h(x2p) - Function_h(x2m)) / 4. / dx1;
    const double m3 = (Function_h(x3p) - Function_h(x3m)) / 6. / dx1;

    return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
}

double MonteCarloEngine::Function_h

(

std::vector< double >

point

)

Definition at line 919 of file MonteCarloEngine.cpp.

                                                           {
    if (point.size() != GetNParameters()) {
        throw std::runtime_error("MCMCENgine::Function_h : Invalid number of entries in the vector.");
    }
    return LogLikelihood(point) + LogAPrioriProbability(point);
}

std::map<std::string, BCH1D * > MonteCarloEngine::getHistograms1D ( ) const

inline

Definition at line 225 of file MonteCarloEngine.h.

    {
        return Histo1D;
    }

std::map<std::string, BCH2D * > MonteCarloEngine::getHistograms2D ( ) const

inline

Definition at line 230 of file MonteCarloEngine.h.

    {
        return Histo2D;
    }

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.

    {
        return HistoLog.str().c_str();
    }

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.

    {
        return NumOfDiscardedEvents;
    }

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.

    {
        return NumOfUsedEvents;
    }

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_i

  the pointer to the model defined in SomeModel.conf

Definition at line 42 of file MonteCarloEngine.cpp.

                                              {
    TH1D * lhisto = new TH1D("LogLikelihood", "LogLikelihood",
                    NBINS1D, 1., -1.);
    lhisto->SetDefaultBufferSize(100000);
    lhisto->GetXaxis()->SetTitle("LogLikelihood");
    BCH1D * bclhisto = new BCH1D(lhisto);
    Histo1D["LogLikelihood"] = bclhisto;
    Mod = Mod_i;
    int k = 0, kweight = 0;
    for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
            it < Obs_ALL.end(); it++) {
        if (!it->isTMCMC()) {
            k++;
            if (it->getDistr().compare("noweight") != 0)
                kweight++;
        }
        std::string HistName = it->getName();
        if (Histo1D.find(HistName) == Histo1D.end()) {
            TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(),
                    NBINS1D, it->getMin(), it->getMax());
            histo->GetXaxis()->SetTitle(it->getLabel().c_str());
            BCH1D * bchisto = new BCH1D(histo);
            Histo1D[HistName] = bchisto;
            thMin[it->getName()] = std::numeric_limits<double>::max();
            thMax[it->getName()] = -std::numeric_limits<double>::max();
        }
    }
    for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
            it < Obs2D_ALL.end(); it++) {
        if ((it->getDistr()).compare("file") == 0) {
            if (!it->isTMCMC())
                throw std::runtime_error("ERROR: cannot handle noMCMC for Observable2D file yet!");
        } else if (it->getDistr().compare("weight") == 0)
            throw std::runtime_error("ERROR: do not use Observable2D for analytic 2D weights!");
        std::string HistName = it->getName();
        if (Histo2D.find(HistName) == Histo2D.end()) {
            TH2D * histo2 = new TH2D(HistName.c_str(),
                    (it->getLabel() + " vs " + it->getLabel2()).c_str(),
                    NBINS2D, it->getMin(), it->getMax(),
                    NBINS2D, it->getMin2(), it->getMax2());
            histo2->GetXaxis()->SetTitle(it->getLabel().c_str());
            histo2->GetYaxis()->SetTitle(it->getLabel2().c_str());
            BCH2D * bchisto2 = new BCH2D(histo2);
            Histo2D[HistName] = bchisto2;
        }
    }
    for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
            it1 != CGO.end(); ++it1) {
        std::vector<Observable> pino(it1->getObs());
        for (std::vector<Observable>::iterator it = pino.begin();
                it != pino.end(); ++it) {
            //for (int i = 0; i < it1->getObs().size(); i++) {
            //Observable * it = &(it1->getObs().at(i));
            if ((it->getDistr()).compare("file") == 0)
                throw std::runtime_error("Cannot use file in CorrelatedGaussianObservables!");
            if (!(it->isTMCMC())) {
                k++;
                if (it->getDistr().compare("noweight") != 0)
                    throw std::runtime_error("Cannot use weight in CorrelatedGaussianObservables!");
            }
            std::string HistName = it->getName();
            if (Histo1D.find(HistName) == Histo1D.end()) {
                TH1D * histo = new TH1D(HistName.c_str(), it->getLabel().c_str(),
                        NBINS1D, it->getMin(), it->getMax());
                histo->GetXaxis()->SetTitle(it->getLabel().c_str());
                BCH1D * bchisto = new BCH1D(histo);
                Histo1D[HistName] = bchisto;
                thMin[HistName] = std::numeric_limits<double>::max();
                thMax[HistName] = -std::numeric_limits<double>::max();
            }
        }
        if (it1->isPrediction()){
            CorrelationMap[it1->getName()] = new TPrincipal(it1->getObs().size());
        }
    }
    kmax = k;
    kwmax = kweight;

    DefineParameters();

};

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]

parameters

the vector of the parameters that are being varied in the Monte Carlo

Returns

the loglikelihood

Definition at line 248 of file MonteCarloEngine.cpp.

                                                                          {
    // This methods returns the logarithm of the conditional probability
    // p(data|parameters). This is where you have to define your model.

    double logprob = 0.;
    
    setDParsFromParameters(parameters, DPars);

    // if update false set probability equal zero
    if (!Mod->Update(DPars)) {
#ifdef _MCDEBUG
        std::cout << "event discarded" << std::endl;

        /* Debug */
        //for (int k = 0; k < parameters.size(); k++)
        //    std::cout << "  " << GetParameter(k)->GetName() << " = "
        //              << DPars[GetParameter(k)->GetName()] << std::endl;
#endif
        NumOfDiscardedEvents++;
        return (log(0.));
    }
#ifdef _MCDEBUG
    //std::cout << "event used in MC" << std::endl;
#endif

    for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it != Obs_ALL.end(); it++) {
        if (it->isTMCMC()) logprob += it->computeWeight();
    }

    for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it != Obs2D_ALL.end(); it++) {
        if (it->isTMCMC()) logprob += it->computeWeight();
    }

    for (std::vector<CorrelatedGaussianObservables>::iterator it = CGO.begin(); it < CGO.end(); it++) {
        if(!(it->isPrediction())) logprob += it->computeWeight();
    }
    if (isnan(logprob)) {
        NumOfDiscardedEvents++;
        std::cout << "Event discarded since logprob evaluated to NAN.\n";
        return (log(0.));
    }
    NumOfUsedEvents++;
    return logprob;
}

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.

                                              {
#ifdef _MPI
    unsigned mychain = 0;
    int iproc = 0;
    unsigned npars = GetNParameters();
    int buffsize = npars + 1;
    int index_chain[procnum];
    double *recvbuff = new double[buffsize];
    std::vector<std::vector<double> > fMCMCxvect;
    std::vector<double> pars;
    double **buff;


    buff = new double*[procnum];
    int obsbuffsize = Obs_ALL.size() + 2 * Obs2D_ALL.size();
    for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++)
        obsbuffsize += it1->getObs().size();
    buff[0] = new double[procnum * obsbuffsize];
    for (int i = 1; i < procnum; i++) {
        buff[i] = buff[i - 1] + obsbuffsize;
        index_chain[i] = -1;
    }

    double ** sendbuff = new double *[procnum];
    sendbuff[0] = new double[procnum * buffsize];
    for (int il = 1; il < procnum; il++)
        sendbuff[il] = sendbuff[il - 1] + buffsize;

    while (mychain < fMCMCNChains) {
        pars.clear();
        for (unsigned int k = 0; k < npars; k++)
            pars.push_back(fMCMCx.at(mychain * npars + k));

        fMCMCxvect.push_back(pars);
        index_chain[iproc] = mychain;
        iproc++;
        mychain++;
        if (iproc < procnum && mychain < fMCMCNChains)
            continue;

        for (int il = 0; il < iproc; il++) {
            //The first entry of the array specifies the task to be executed.

            sendbuff[il][0] = 2.; // 2 = observables calculation
            for (int im = 1; im < buffsize; im++)
                sendbuff[il][im] = fMCMCxvect[il][im - 1];
        }
        for (int il = iproc; il < procnum; il++) {
            sendbuff[il][0] = 3.; // 3 = nothing to execute, but return a buffer of observables
            index_chain[il] = -1;
        }
        //       double inittime = MPI::Wtime();
        MPI::COMM_WORLD.Scatter(sendbuff[0], buffsize, MPI::DOUBLE,
                recvbuff, buffsize, MPI::DOUBLE,
                0);

        if (recvbuff[0] == 2.) { // compute observables
            double sbuff[obsbuffsize];
            std::map<std::string, double> DPars;
            pars.assign(recvbuff + 1, recvbuff + buffsize);
            setDParsFromParameters(pars,DPars);
            Mod->Update(DPars);
          
            int k = 0;
            for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end(); it++) {
                sbuff[k++] = it->computeTheoryValue();
            }
            for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin(); it < Obs2D_ALL.end(); it++) {
                sbuff[k++] = it->computeTheoryValue();
                sbuff[k++] = it->computeTheoryValue2();
            }

            for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin(); it1 < CGO.end(); it1++) {
                std::vector<Observable> pino(it1->getObs());
                for (std::vector<Observable>::iterator it = pino.begin(); it != pino.end(); ++it)
                    sbuff[k++] = it->computeTheoryValue();
            }
            MPI::COMM_WORLD.Gather(sbuff, obsbuffsize, MPI::DOUBLE, buff[0], obsbuffsize, MPI::DOUBLE, 0);
        } else if (recvbuff[0] == 3.) { // do not compute observables, but gather the buffer
            double sbuff[obsbuffsize];
            MPI::COMM_WORLD.Gather(sbuff, obsbuffsize, MPI::DOUBLE, buff[0], obsbuffsize, MPI::DOUBLE, 0);
        }

        for (int il = 0; il < procnum; il++) {
            if (index_chain[il] >= 0) {
                int k = 0;
                // fill the histograms for observables
                int ko = 0, kweight = 0;
                for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
                        it < Obs_ALL.end(); it++) {
                    double th = buff[il][k++];
                    /* set the min and max of theory values */
                    if (th < thMin[it->getName()]) thMin[it->getName()] = th;
                    if (th > thMax[it->getName()]) thMax[it->getName()] = th;
                    Histo1D[it->getName()]->GetHistogram()->Fill(th);
                    if (!it->isTMCMC()) {
                        obval[index_chain[il] * kmax + ko++] = th;
                        if (it->getDistr().compare("noweight") != 0) {
                            obweight[index_chain[il] * kwmax + kweight++] = it->computeWeight(th);
                        }
                    }
                }

                // fill the 2D histograms for observables
                for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
                        it < Obs2D_ALL.end(); it++) {
                    double th1 = buff[il][k++];
                    double th2 = buff[il][k++];
                    Histo2D[it->getName()]->GetHistogram()->Fill(th1, th2);
                }

                // fill the histograms for correlated observables
                for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
                        it1 < CGO.end(); it1++) {
                    std::vector<Observable> pino(it1->getObs());
                    Double_t * COdata = new Double_t[pino.size()];
                    int nObs = 0;
                    for (std::vector<Observable>::iterator it = pino.begin();
                            it != pino.end(); ++it) {
                        double th = buff[il][k++];
                        /* set the min and max of theory values */
                        if (th < thMin[it->getName()]) thMin[it->getName()] = th;
                        if (th > thMax[it->getName()]) thMax[it->getName()] = th;
                        Histo1D[it->getName()]->GetHistogram()->Fill(th);
                        if (it1->isPrediction()) COdata[nObs++] = th;
                    }
                    if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
                    delete [] COdata;
                }
            }
        }
        iproc = 0;
        fMCMCxvect.clear();
    }
    delete sendbuff[0];
    delete [] sendbuff;
    delete [] recvbuff;
    delete buff[0];
    delete [] buff;
#else
    for (unsigned int i = 0; i < fMCMCNChains; ++i) {
        std::vector<double>::const_iterator first = fMCMCx.begin() + i * GetNParameters();
        std::vector<double>::const_iterator last = first + GetNParameters();
        std::vector<double> currvec(first, last);
        setDParsFromParameters(currvec,DPars);

        Mod->Update(DPars);

        // fill the histograms for observables
        int k = 0, kweight = 0;
        for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin();
                it < Obs_ALL.end(); it++) {
            double th = it->computeTheoryValue();
            /* set the min and max of theory values */
            if (th < thMin[it->getName()]) thMin[it->getName()] = th;
            if (th > thMax[it->getName()]) thMax[it->getName()] = th;
            Histo1D[it->getName()]->GetHistogram()->Fill(th);
            if (!it->isTMCMC()) {
                obval[i * kmax + k] = th;
                k++;
                if (it->getDistr().compare("noweight") != 0) {
                    obweight[i * kwmax + kweight] = it->computeWeight();
                    kweight++;
                }
            }
        }

        // fill the 2D histograms for observables
        for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
                it < Obs2D_ALL.end(); it++) {
            double th1 = it->computeTheoryValue();
            double th2 = it->computeTheoryValue2();
            Histo2D[it->getName()]->GetHistogram()->Fill(th1, th2);
        }

        // fill the histograms for correlated observables
        for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
                it1 < CGO.end(); it1++) {
            std::vector<Observable> pino(it1->getObs());
            Double_t * COdata = new Double_t[pino.size()];
            int nObs = 0;
            for (std::vector<Observable>::iterator it = pino.begin();
                    it != pino.end(); ++it) {
                double th = it->computeTheoryValue();
                /* set the min and max of theory values */
                if (th < thMin[it->getName()]) thMin[it->getName()] = th;
                if (th > thMax[it->getName()]) thMax[it->getName()] = th;
                Histo1D[it->getName()]->GetHistogram()->Fill(th);
                if (it1->isPrediction()) COdata[nObs++] = th;
            }
            if (it1->isPrediction()) CorrelationMap[it1->getName()]->AddRow(COdata);
            delete [] COdata;
        }
    }
#endif
    for (unsigned int i = 0; i < fMCMCNChains; i++)
    {
        Histo1D["LogLikelihood"]->GetHistogram()->Fill(MCMCGetLogProbx(i)-LogAPrioriProbability(MCMCGetx(i)));
    }
}

void MonteCarloEngine::PrintCorrelationMatrix

(

const std::string

filename

)

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

Parameters

[in]

filename

the name of the file where the correlation matrix is printed

Definition at line 611 of file MonteCarloEngine.cpp.

                                                                      {
    std::ofstream out;
    out.open(filename.c_str(), std::ios::out);

    int npar = GetNParameters();

    for (int i = 0; i < npar; ++i)
        out << " & " << GetParameter(i)->GetName();
    out << " \\\\" << std::endl;

    for (int i = 0; i < npar; ++i) {
        out << GetParameter(i)->GetName() << " & $";
        for (int j = 0; j < npar; ++j) {
            if (i != j) {
                BCH2D* bch2d_temp = GetMarginalized(GetParameter(i), GetParameter(j));
                if (bch2d_temp != NULL)
                    out << bch2d_temp->GetHistogram()->GetCorrelationFactor();
                else
                    out << 0.;
            } else
                out << 1.;
            if (j == npar - 1) out << "$ \\\\" << std::endl;
            else out << "$ & $";
        }
    }

    out.close();
}

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]	out	a reference to an object of type BCModelOutput as defined in the BAT libraries
[in]	it	a reference to an object of type Observable
[in]	OutputDir	the name of the output directory

Definition at line 520 of file MonteCarloEngine.cpp.

                                   {
    std::string HistName = it.getName();
    double min = thMin[it.getName()];
    double max = thMax[it.getName()];

    if (Histo1D[HistName]->GetHistogram()->Integral() > 0.0) {
        std::string fname = OutputDir + "/" + HistName + ".pdf";
        //        BCH1D* pippo =  Histo1D[HistName];
        //        double x = pippo->GetMean();
        //        pippo->Print("Dmd1.pdf");
        Histo1D[HistName]->SetGlobalMode(it.computeTheoryValue());
        Histo1D[HistName]->Print(fname.c_str());
        std::cout << fname << " has been created." << std::endl;
        out.Write(Histo1D[HistName]->GetHistogram());
        CheckHistogram(*Histo1D[HistName]->GetHistogram(), it.getName());
    } else
        HistoLog << "WARNING: The histogram of "
            << it.getName() << " is empty!" << std::endl;

    HistoLog.precision(10);
    HistoLog << "  [min, max]=[" << min << ", " << max << "]" << std::endl;
    HistoLog.precision(6);
}

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]	out	a reference to an object of type BCModelOutput as defined in the BAT libraries
[in]	OutputDir	the name of the output directory

Definition at line 545 of file MonteCarloEngine.cpp.

                                                                                    {
    std::vector<double> mode(GetBestFitParameters());
    if (mode.size() == 0) {
        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");
        else sleep(2);
    }
    setDParsFromParameters(mode,DPars);

    Mod->Update(DPars);

    // print the histograms to pdf files
    for (boost::ptr_vector<Observable>::iterator it = Obs_ALL.begin(); it < Obs_ALL.end();
            it++) {
        PrintHistogram(out, *it, OutputDir);
    }
    for (std::vector<CorrelatedGaussianObservables>::iterator it1 = CGO.begin();
            it1 < CGO.end(); it1++) {
        std::vector<Observable> pino(it1->getObs());
        for (std::vector<Observable>::iterator it = pino.begin();
                it != pino.end(); ++it)
            PrintHistogram(out, *it, OutputDir);
    }
    for (std::vector<Observable2D>::iterator it = Obs2D_ALL.begin();
            it < Obs2D_ALL.end(); it++) {
        std::string HistName = it->getName();
        if (Histo2D[HistName]->GetHistogram()->Integral() > 0.0) {
            std::string fname = OutputDir + "/" + HistName + ".pdf";
            double th[2];
            th[0] = it->computeTheoryValue();
            th[1] = it->computeTheoryValue2();
            Histo2D[HistName]->SetGlobalMode(th);
            Histo2D[HistName]->Print(fname.c_str());
            std::cout << fname << " has been created." << std::endl;
            out.Write(Histo2D[HistName]->GetHistogram());
            CheckHistogram(*Histo2D[HistName]->GetHistogram(), HistName);
        } else
            HistoLog << "WARNING: The histogram of "
                << HistName << " is empty!" << std::endl;
    }
    std::string fname = OutputDir + "/LogLikelihood.pdf";
    Histo1D["LogLikelihood"]->Print(fname.c_str());
    std::cout << fname << " has been created." << std::endl;
    out.Write(Histo1D["LogLikelihood"]->GetHistogram());
}

double MonteCarloEngine::SecondDerivative	(	const BCParameter *	par1,
		const BCParameter *	par2,
		std::vector< double >	point
	)

Definition at line 849 of file MonteCarloEngine.cpp.

                                                                                                                     {

    if (point.size() != GetNParameters()) {
        throw std::runtime_error("MCMCENgine::SecondDerivative : Invalid number of entries in the vector.");
    }

    // define steps
    const double dy1 = par2->GetRangeWidth() / NSTEPS;
    const double dy2 = dy1 * 2.;
    const double dy3 = dy1 * 3.;

    // define points at which to evaluate
    std::vector<double> y1p = point;
    std::vector<double> y1m = point;
    std::vector<double> y2p = point;
    std::vector<double> y2m = point;
    std::vector<double> y3p = point;
    std::vector<double> y3m = point;

    unsigned idy = GetParameters().Index(par2->GetName());

    y1p[idy] += dy1;
    y1m[idy] -= dy1;
    y2p[idy] += dy2;
    y2m[idy] -= dy2;
    y3p[idy] += dy3;
    y3m[idy] -= dy3;

    const double m1 = (FirstDerivative(par1, y1p) - FirstDerivative(par1, y1m)) / 2. / dy1;
    const double m2 = (FirstDerivative(par1, y2p) - FirstDerivative(par1, y2m)) / 4. / dy1;
    const double m3 = (FirstDerivative(par1, y3p) - FirstDerivative(par1, y3m)) / 6. / dy1;

    return 3. / 2. * m1 - 3. / 5. * m2 + 1. / 10. * m3;
}

void MonteCarloEngine::setDParsFromParameters	(	const std::vector< double > &	parameters,
		std::map< std::string, double > &	DPars_i
	)

Definition at line 199 of file MonteCarloEngine.cpp.

{
    std::map<std::string, std::vector<double> > cgpmap;

    unsigned int k = 0;
    for (std::vector<ModelParameter>::const_iterator it = ModPars.begin(); it < ModPars.end(); it++){
        if(it->IsFixed())
            continue;
        if(it->getname().compare(GetParameter(k)->GetName()) != 0)
            {
                        std::stringstream out;
                        out << it->getname();
                        throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + "is sitting at the wrong position in the BAT parameters vector");
                    }
        if (it->IsCorrelated()) {
            std::string index = it->getname().substr(it->getCgp_name().size());
            unsigned long int lindex = strtol(index.c_str(),NULL,10);
            if (lindex - 1 == cgpmap[it->getCgp_name()].size())
                cgpmap[it->getCgp_name()].push_back(parameters[k]);
            else {
                std::stringstream out;
                out << it->getname() << " " << lindex;
                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");
            }

        } else
            DPars_i[it->getname()] = parameters[k];
        k++;
    }

    for (unsigned int j = 0; j < CGP.size(); j++) {
        std::vector<double> current = cgpmap.at(CGP[j].getName());
        if (current.size() != CGP[j].getPars().size()) {
            std::stringstream out;
            out << CGP[j].getName();
            throw std::runtime_error("MonteCarloEngine::setDParsFromParameters(): " + out.str() + " appears to be represented in cgpmap with a wrong size");
        }
        
        std::vector<double> porig = CGP[j].getOrigParsValue(current);

        for(unsigned int l = 0; l < porig.size(); l++) {
            DPars_i[CGP[j].getPar(l).getname()] = porig[l];
        }
    }
}

void MonteCarloEngine::setHistogramOverFlow ( bool  overflow )

inline

Definition at line 235 of file MonteCarloEngine.h.

    {
        if (overflow) TH1::StatOverflows(kTRUE);
        else TH1::StatOverflows(kFALSE);
    }

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]

i

the number of chains

Definition at line 124 of file MonteCarloEngine.cpp.

                                                {
    MCMCSetNChains(i);
    obval = new double[fMCMCNChains * kmax];
    obweight = new double[fMCMCNChains * kwmax];
}

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.

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

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:

• MonteCarloEngine.h
• MonteCarloEngine.cpp