BCEngineMCMC.h

#ifndef __BCENGINEMCMC__H
#define __BCENGINEMCMC__H

/*
 * Copyright (C) 2007-2013, the BAT core developer team
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 * For documentation see http://mpp.mpg.de/bat
 */

// ---------------------------------------------------------

#include "BCH1D.h"
#include "BCH2D.h"
#include "BCLog.h"
#include "BCObservable.h"
#include "BCParameter.h"
#include "BCParameterSet.h"

#include <TMatrixD.h>
#include <TMatrixDSym.h>
#include <TRandom3.h>
#include <TVectorD.h>

#include <algorithm>
#include <limits>
#include <string>
#include <utility>
#include <vector>

// ROOT classes
class TF1;
class TFile;
class TH1;
class TH2;
class TTree;

class BCVariable;

// ---------------------------------------------------------

class BCEngineMCMC
{

public:

    enum Precision {
        kLow,
        kQuick,
        kMedium,
        kHigh,
        kVeryHigh
    };

    enum Phase {
        kPreRun     = -1, 
        kUnsetPhase =  0, 
        kMainRun    = +1  
    };

    enum InitialPositionScheme {
        kInitCenter           =  0, 
        kInitRandomUniform    =  1, 
        kInitUserDefined      =  2, 
        kInitRandomPrior      =  3  
    };

private:
    // grant friendship to BCModel to access the ThreadLocalStorage
    friend class BCModel;

    struct ThreadLocalStorage {
        std::vector<double> parameters;

        double log_prior;

        double log_likelihood;

        double log_probability;

        TRandom3* rng;

        TVectorD yLocal;

        ThreadLocalStorage(unsigned dim);

        ThreadLocalStorage(const ThreadLocalStorage& other);

        ThreadLocalStorage& operator = (ThreadLocalStorage other);

        void swap(BCEngineMCMC::ThreadLocalStorage& A, BCEngineMCMC::ThreadLocalStorage& B);

        virtual ~ThreadLocalStorage();

        double scale(double dof);
    };

public:

    struct ChainState {
        unsigned iteration;     
        std::vector<double> parameters; 
        std::vector<double> observables; 
        double log_probability; 
        double log_likelihood;  
        double log_prior;       

        ChainState(int n_obs = 0) : iteration(0),
            observables(std::vector<double>(n_obs, 0.)),
            log_probability(-std::numeric_limits<double>::infinity()),
            log_likelihood(-std::numeric_limits<double>::infinity()),
            log_prior(-std::numeric_limits<double>::infinity())
        {}

        ChainState& operator=(const ThreadLocalStorage& tls)
        {
            parameters = tls.parameters;
            log_prior = tls.log_prior;
            log_likelihood = tls.log_likelihood;
            log_probability = tls.log_probability;
            return *this;
        }
    };

    struct Statistics {

        Statistics(unsigned n_par = 0, unsigned n_obs = 0);

        unsigned n_samples;                           
        std::vector<double> mean;                     
        std::vector<double> variance;                 
        std::vector<double> stderrpar;                
        std::vector<double> stderrobs;                
        std::vector<std::vector<double> > covariance; 
        std::vector<double> minimum;                  
        std::vector<double> maximum;                  
        double probability_mean;                      
        double probability_variance;                  
        std::vector<double> modepar;                  
        std::vector<double> modeobs;                  
        double probability_at_mode;                   
        unsigned n_samples_efficiency;                
        std::vector<double> efficiency;               

        void Clear(bool clear_mode = true, bool clear_efficiency = true);

        void Init(unsigned n_par, unsigned n_obs);

        void Reset(bool reset_mode = true, bool reset_efficiency = true);

        void ResetEfficiencies();

        Statistics& operator += (const Statistics& rhs);

        void Update(const ChainState& cs);
    };

    BCEngineMCMC(const std::string& name = "model");

    BCEngineMCMC(const BCEngineMCMC& enginemcmc);

    BCEngineMCMC(const std::string& filename, const std::string& name, bool loadObservables = true);

    BCEngineMCMC& operator=(const BCEngineMCMC&);

    virtual ~BCEngineMCMC();

    const std::string& GetName() const
    { return fName; }

    const std::string& GetSafeName() const
    { return fSafeName; }

    unsigned GetNChains() const
    { return fMCMCNChains; }

    unsigned GetNLag() const
    { return fMCMCNLag; }

    int GetCurrentIteration() const
    { return fMCMCCurrentIteration; }

    unsigned GetCurrentChain() const;

    int GetNIterationsConvergenceGlobal() const
    { return fMCMCNIterationsConvergenceGlobal; }

    unsigned GetNIterationsPreRun() const;

    unsigned GetNIterationsPreRunMin() const
    { return fMCMCNIterationsPreRunMin; }

    unsigned GetNIterationsPreRunMax() const
    { return fMCMCNIterationsPreRunMax; }

    unsigned GetNIterationsRun() const
    { return fMCMCNIterationsRun; }

    unsigned GetNIterationsPreRunCheck() const
    { return fMCMCNIterationsPreRunCheck; }

    unsigned GetPreRunCheckClear()
    { return fMCMCPreRunCheckClear; }

    double GetMinimumEfficiency() const
    { return fMCMCEfficiencyMin; }

    double GetMaximumEfficiency() const
    { return fMCMCEfficiencyMax; }

    double GetScaleFactorLowerLimit() const
    { return fMCMCScaleFactorLowerLimit; }

    double GetScaleFactorUpperLimit() const
    { return fMCMCScaleFactorUpperLimit; }

    const std::vector<std::vector<double> >& GetScaleFactors() const
    { return fMCMCProposalFunctionScaleFactor; }

    const ChainState& GetChainState(unsigned c) const
    { return fMCMCStates.at(c); }

    const std::vector<double>& Getx(unsigned c) const
    { return fMCMCStates.at(c).parameters; }

    double Getx(unsigned c, unsigned p) const
    { return fMCMCStates.at(c).parameters.at(p); }

    double GetLogProbx(unsigned c) const
    { return fMCMCStates.at(c).log_probability; }

    BCEngineMCMC::Phase GetPhase() const
    { return fMCMCPhase; }

    BCEngineMCMC::InitialPositionScheme GetInitialPositionScheme() const
    { return fInitialPositionScheme; }

    unsigned GetInitialPositionAttemptLimit() const
    { return fInitialPositionAttemptLimit; }

    bool GetProposeMultivariate() const
    { return fMCMCProposeMultivariate; }

    double GetProposalFunctionDof() const
    { return fMCMCProposalFunctionDof; }

    unsigned GetMultivariateCovarianceUpdates() const
    { return fMultivariateCovarianceUpdates; }

    /*
     * @return weighting parameter for multivariate proposal function covariance update. */
    double GetMultivariateCovarianceUpdateLambda() const
    { return fMultivariateCovarianceUpdateLambda; }

    double GetMultivariateEpsilon() const
    { return fMultivariateEpsilon; }

    double GetMultivariateScaleMultiplier() const
    { return fMultivariateScaleMultiplier; }

    double GetRValueParametersCriterion() const
    { return fMCMCRValueParametersCriterion; }

    const std::vector<double>& GetRValueParameters() const
    { return fMCMCRValueParameters; }

    double GetRValueParameters(unsigned index) const
    { return fMCMCRValueParameters.at(index); }

    bool GetCorrectRValueForSamplingVariability() const
    { return fCorrectRValueForSamplingVariability; }

    bool GetFlagRun() const
    { return fMCMCPhase == kMainRun; }

    TTree* GetMarkovChainTree() const
    { return fMCMCTree;}

    TTree* GetParameterTree() const
    { return fParameterTree;}

    TFile* GetOutputFile() const
    { return fMCMCOutputFile; }

    const BCEngineMCMC::Statistics& GetStatistics() const
    { return fMCMCStatistics_AllChains; }

    const BCEngineMCMC::Statistics& GetStatistics(unsigned c) const
    { return fMCMCStatistics.at(c); }

    const std::vector<BCEngineMCMC::Statistics>& GetStatisticsVector() const
    { return fMCMCStatistics; }

    bool GetRescaleHistogramRangesAfterPreRun() const
    { return fRescaleHistogramRangesAfterPreRun; }

    double GetHistogramRescalePadding() const
    { return fHistogramRescalePadding; }

    virtual std::vector<unsigned> GetH1DPrintOrder() const;

    virtual std::vector<std::pair<unsigned, unsigned> > GetH2DPrintOrder() const;

    bool MarginalizedHistogramExists(unsigned index) const
    { return index < fH1Marginalized.size() and fH1Marginalized[index]; }

    bool MarginalizedHistogramExists(unsigned index1, unsigned index2) const
    { return index1 < fH2Marginalized.size() and index2 < fH2Marginalized[index1].size() and fH2Marginalized[index1][index2]; }

    TH1* GetMarginalizedHistogram(const std::string& name) const
    { return GetMarginalizedHistogram(fParameters.Index(name)); }

    TH1* GetMarginalizedHistogram(unsigned index) const;

    TH2* GetMarginalizedHistogram(const std::string& name1, const std::string& name2) const
    { return GetMarginalizedHistogram(fParameters.Index(name1), fParameters.Index(name2)); }

    TH2* GetMarginalizedHistogram(unsigned index1, unsigned index2) const;

    BCH1D GetMarginalized(const std::string& name) const
    { return GetMarginalized(fParameters.Index(name)); }

    BCH1D GetMarginalized(unsigned index) const;

    BCH2D GetMarginalized(const std::string& name1, const std::string& name2) const
    { return GetMarginalized(fParameters.Index(name1), fParameters.Index(name2)); }

    BCH2D GetMarginalized(unsigned index1, unsigned index2) const;

    unsigned GetMaximumParameterNameLength(bool observables = true) const
    { return (observables) ? std::max(fParameters.MaxNameLength(), fObservables.MaxNameLength()) : fParameters.MaxNameLength(); }

    BCVariable& GetVariable(unsigned index)
    { return (index < GetNParameters()) ? static_cast<BCVariable&>(fParameters[index]) : static_cast<BCVariable&>(fObservables.At(index - GetNParameters())); }

    const BCVariable& GetVariable(unsigned index) const
    { return (index < GetNParameters()) ? static_cast<const BCVariable&>(fParameters[index]) : static_cast<const BCVariable&>(fObservables.At(index - GetNParameters())); }

    unsigned GetNVariables() const
    { return fParameters.Size() + fObservables.Size(); }

    BCParameterSet& GetParameters()
    { return fParameters; }

    const BCParameterSet& GetParameters() const
    { return fParameters; }

    BCParameter& GetParameter(unsigned index)
    { return fParameters.At(index); }

    const BCParameter& GetParameter(unsigned index) const
    { return fParameters.At(index); }

    BCParameter& GetParameter(const std::string& name)
    { return fParameters.Get(name); }

    const BCParameter& GetParameter(const std::string& name) const
    { return fParameters.Get(name); }

    unsigned GetNParameters() const
    { return fParameters.Size(); }

    unsigned GetNFixedParameters() const
    { return fParameters.GetNFixedParameters(); }

    unsigned GetNFreeParameters() const
    { return fParameters.GetNFreeParameters(); }

    BCObservableSet& GetObservables()
    { return fObservables; }

    const BCObservableSet& GetObservables() const
    { return fObservables; }

    BCObservable& GetObservable(unsigned index)
    { return fObservables.At(index); }

    const BCObservable& GetObservable(unsigned index) const
    { return fObservables.At(index); }

    BCObservable& GetObservable(const std::string& name)
    { return fObservables.Get(name); }

    const BCObservable& GetObservable(const std::string& name) const
    { return fObservables.Get(name); }

    unsigned GetNObservables() const
    { return fObservables.Size(); }

    virtual const std::vector<double>& GetBestFitParameters() const;

    virtual const std::vector<double>& GetBestFitParameterErrors() const;

    const std::vector<double>& GetLocalModes(bool force_recalculation = false);

    virtual double GetLogMaximum() const
    { return fMCMCStatistics_AllChains.probability_at_mode; };

    bool GetReuseObservables() const
    { return fMCMCTreeReuseObservables; }

    BCH1D& GetBCH1DdrawingOptions()
    { return fBCH1DdrawingOptions; }

    BCH2D& GetBCH2DdrawingOptions()
    { return fBCH2DdrawingOptions; }

    void SetName(const std::string& name);

    void SetScaleFactorLowerLimit(double l)
    { fMCMCScaleFactorLowerLimit = l; }

    void SetScaleFactorUpperLimit(double l)
    { fMCMCScaleFactorUpperLimit = l; }

    void SetInitialScaleFactors(const std::vector<double>& scale)
    { fMCMCInitialScaleFactors = scale; }

    void SetNChains(unsigned n)
    { fMCMCNChains = n; }

    void SetNLag(unsigned n);

    void SetNIterationsPreRunMax(unsigned n)
    { fMCMCNIterationsPreRunMax = n; }

    void SetNIterationsRun(unsigned n)
    { fMCMCNIterationsRun = n; }

    void SetNIterationsPreRunMin(unsigned n)
    { fMCMCNIterationsPreRunMin = n; }

    void SetNIterationsPreRunCheck(unsigned n)
    { fMCMCNIterationsPreRunCheck = n; }

    void SetPreRunCheckClear(unsigned n)
    { fMCMCPreRunCheckClear = n; }

    void SetMinimumEfficiency(double efficiency)
    { fMCMCEfficiencyMin = efficiency; }

    void SetMaximumEfficiency(double efficiency)
    { fMCMCEfficiencyMax = efficiency; }

    void SetRandomSeed(unsigned seed);

    void SetInitialPositions(const std::vector<double>& x0s);

    void SetInitialPositions(const std::vector< std::vector<double> >& x0s)
    { fMCMCInitialPosition = x0s; SetInitialPositionScheme(BCEngineMCMC::kInitUserDefined); }

    void SetInitialPositionScheme(BCEngineMCMC::InitialPositionScheme scheme)
    { fInitialPositionScheme = scheme; }

    void SetInitialPositionAttemptLimit(unsigned n)
    { fInitialPositionAttemptLimit = n; }

    void SetProposeMultivariate(bool flag)
    {
        if (flag != fMCMCProposeMultivariate && fMCMCPhase != kUnsetPhase)
            BCLog::OutWarning("You have changed the proposal function flag, but not required the prerun be rerun.");
        fMCMCProposeMultivariate = flag;
    }

    void SetProposalFunctionDof(double dof = 1)
    {
        fMCMCProposalFunctionDof = dof;
    }

    void SetMultivariateCovarianceUpdateLambda(double l)
    { fMultivariateCovarianceUpdateLambda = std::max<double>(0, std::min<double>(1, l)); }

    void SetMultivariateEpsilon(double epsilon)
    { fMultivariateEpsilon = epsilon; }

    void SetMultivariateScaleMultiplier(double s)
    { fMultivariateScaleMultiplier = s; }

    void SetFlagFillHistograms(bool flag)
    { fParameters.FillHistograms(flag); fObservables.FillHistograms(flag); }

    void SetFlagFillHistograms(bool flag_1d, bool flag_2d)
    { fParameters.FillHistograms(flag_1d, flag_2d); fObservables.FillHistograms(flag_1d, flag_2d); }

    void SetFillHistogramParPar(unsigned x, unsigned y, bool flag = true)
    { SetFillHistogram(x, y, flag); }

    void SetFillHistogramParPar(const std::string& x, const std::string& y, bool flag = true)
    { SetFillHistogramParPar(fParameters.Index(x), fParameters.Index(y), flag); }

    void SetFillHistogramParObs(unsigned x, unsigned y, bool flag = true)
    { SetFillHistogram(x, -(y + 1), flag); }

    void SetFillHistogramParObs(const std::string& x, const std::string& y, bool flag = true)
    { SetFillHistogramParObs(fParameters.Index(x), fObservables.Index(y), flag); }

    void SetFillHistogramObsObs(unsigned x, unsigned y, bool flag = true)
    { SetFillHistogram(-(x + 1), -(y + 1), flag); }

    void SetFillHistogramObsObs(const std::string& x, const std::string& y, bool flag = true)
    { SetFillHistogramObsObs(fObservables.Index(x), fObservables.Index(y), flag); }

    void SetFillHistogramObsPar(unsigned x, unsigned y, bool flag = true)
    { SetFillHistogram(-(x + 1), y, flag); }

    void SetFillHistogramObsPar(const std::string& x, const std::string& y, bool flag = true)
    { SetFillHistogramObsPar(fObservables.Index(x), fParameters.Index(y), flag); }

    void SetFlagPreRun(bool flag)
    { if (flag) fMCMCPhase = kUnsetPhase; }

    void SetRValueParametersCriterion(double r)
    { fMCMCRValueParametersCriterion = r; }

    void SetCorrectRValueForSamplingVariability(bool flag = true)
    { fCorrectRValueForSamplingVariability = flag; }

    void SetPrecision(BCEngineMCMC::Precision precision);

    void SetPrecision(const BCEngineMCMC* other)
    { if (other) SetPrecision(*other); }

    void SetPrecision(const BCEngineMCMC& other);

    void SetNbins(unsigned int nbins)
    { fParameters.SetNBins(nbins); fObservables.SetNBins(nbins); }

    void SetReuseObservables(bool flag)
    { fMCMCTreeReuseObservables = flag; }

    void SetRescaleHistogramRangesAfterPreRun(bool flag = true)
    { fRescaleHistogramRangesAfterPreRun = flag; }

    void SetHistogramRescalingPadding(double factor)
    { fHistogramRescalePadding = factor; }

    void WriteMarkovChain(bool flag)
    { WriteMarkovChainRun(flag); WriteMarkovChainPreRun(flag); }

    void WriteMarkovChainRun(bool flag);

    void WriteMarkovChainPreRun(bool flag);

    void WriteMarkovChain(const std::string& filename, const std::string& option, bool flag_run = true, bool flag_prerun = true);

    void SetPriorConstant(unsigned index)
    { fParameters.At(index).SetPriorConstant(); }

    void SetPriorConstant(const std::string& name)
    { SetPriorConstant(fParameters.Index(name)); }

    void SetPrior(unsigned index, TF1& f, bool logL = true);

    void SetPrior(const std::string& name, TF1& f, bool logL = true)
    { return SetPrior(fParameters.Index(name), f, logL); }

    void SetPriorDelta(unsigned index, double value)
    { fParameters.At(index).Fix(value); }

    void SetPriorDelta(const std::string& name, double value)
    { SetPriorDelta(fParameters.Index(name), value); }

    void SetPriorGauss(unsigned index, double mean, double sigma);

    void SetPriorGauss(const std::string& name, double mean, double sigma)
    { SetPriorGauss(fParameters.Index(name), mean, sigma); }

    void SetPriorGauss(unsigned index, double mode, double sigma_below, double sigma_above);

    void SetPriorGauss(const std::string& name, double mode, double sigma_below, double sigma_above)
    {   SetPriorGauss(fParameters.Index(name), mode, sigma_below, sigma_above); }

    void SetPrior(unsigned index, TH1& h, bool interpolate = false);

    void SetPrior(const std::string& name, TH1& h, bool interpolate = false)
    { SetPrior(fParameters.Index(name), h, interpolate); }

    void SetPriorConstantAll()
    { fParameters.SetPriorConstantAll(); }

    void WriteMarginalizedDistributions(const std::string& filename, const std::string& option, bool closeExistingFile = false);

    virtual void PrintSummary() const;

    void PrintParameters(const std::vector<double>& P, void (*output)(const std::string&) = BCLog::OutSummary) const;

    unsigned PrintAllMarginalized(const std::string& filename, unsigned hdiv = 1, unsigned vdiv = 1) const;

    unsigned PrintParameterPlot(const std::string& filename, int npar = 10, double interval_content = 68e-2, std::vector<double> quantile_values = std::vector<double>(0), bool rescale_ranges = true) const;

    bool DrawParameterPlot(unsigned i0, unsigned npar = 0, double interval_content = 68e-2, std::vector<double> quantile_values = std::vector<double>(0), bool rescale_ranges = true) const;

    bool PrintCorrelationMatrix(const std::string& filename = "matrix.pdf") const;

    bool PrintCorrelationPlot(const std::string& filename = "correlation.pdf", bool include_observables = true) const;

    bool PrintParameterLatex(const std::string& filename) const;

    virtual void CreateHistograms(bool rescale_ranges = false);

    virtual void InitializeMarkovChainTree(bool replacetree = false, bool replacefile = false);

    virtual bool AddParameter(const std::string& name, double min, double max, const std::string& latexname = "", const std::string& unitstring = "")
    { return fParameters.Add(name, min, max, latexname, unitstring); }

    virtual bool AddParameter(BCParameter& parameter)
    { return fParameters.Add(parameter); }

    virtual bool AddObservable(const std::string& name, double min, double max, const std::string& latexname = "", const std::string& unitstring = "")
    { return fObservables.Add(name, min, max, latexname, unitstring); }

    virtual bool AddObservable(BCObservable& obs)
    { return fObservables.Add(obs); }

    virtual void EvaluateObservables();

    virtual void EvaluateObservables(unsigned chain);

    virtual void CalculateObservables(const std::vector<double>& pars)
    { (void)pars; }

    virtual double ProposalFunction(unsigned ichain, unsigned ipar);

    bool GetProposalPointMetropolis(unsigned chain, std::vector<double>& x);

    bool GetProposalPointMetropolis(unsigned chain, unsigned parameter, std::vector<double>& x);

    bool GetNewPointMetropolis();

    bool GetNewPointMetropolis(unsigned chain);

    bool GetNewPointMetropolis(unsigned chain, unsigned parameter);

    bool AcceptOrRejectPoint(unsigned chain, unsigned parameter);

    void InChainFillHistograms(const ChainState& cs);

    void InChainFillHistograms();

    void InChainFillTree(const ChainState& cs, unsigned chain_number);

    void InChainFillTree();

    virtual double LogEval(const std::vector<double>& parameters) = 0;

    bool Metropolis();

    bool MetropolisPreRun();

    static double RValue(const std::vector<double>& means, const std::vector<double>& variances, unsigned n, bool correctForSamplingVariability = true);

    void MCMCInitialize();

    virtual void MCMCUserInitialize()
    {}

    virtual void ResetResults();

    virtual void MCMCUserIterationInterface()
    {}

    virtual void MCMCCurrentPointInterface(const std::vector<double>& point, int ichain, bool accepted)
    {
        (void)point;
        (void)ichain;
        (void)accepted;
    }

    void LoadParametersFromTree(TTree* partree, bool loadObservables = true);

    void LoadMCMCParameters(TTree& partree);

    virtual bool ParameterTreeMatchesModel(TTree* partree, bool checkObservables = true);

    void LoadMCMC(const std::string& filename, std::string mcmcTreeName = "", std::string parameterTreeName = "", bool loadObservables = true);

    void LoadMCMC(TTree* mcmcTree, TTree* parTree, bool loadObservables = true);

    bool ValidMCMCTree(TTree* tree, bool checkObservables = true) const;

    bool ValidParameterTree(TTree* tree) const;

    void CloseOutputFile();

    virtual void Remarginalize(bool autorange = true);

    void PrepareToContinueMarginalization(const std::string& filename, const std::string& mcmcTreeName = "", const std::string& parameterTreeName = "", bool loadObservables = true, bool autorange = true);

    virtual bool UpdateMultivariateProposalFunctionCovariances(double a);

    virtual bool UpdateMultivariateProposalFunctionCovariances();

    void CalculateCholeskyDecompositions();

    void UpdateChainIndex(int chain);

private:

    std::vector<ThreadLocalStorage> fMCMCThreadLocalStorage;

    void SyncThreadStorage();

    typedef std::vector<unsigned> ChainIndex_t;
    ChainIndex_t fChainIndex;

    void DeleteMarginals();

    void CloneMarginals(const BCEngineMCMC&);

protected:

    void SetFillHistogram(int x, int y, bool flag);

    virtual void PrintModelSummary() const;

    virtual void PrintBestFitSummary() const;

    virtual std::string GetBestFitSummary(unsigned i) const;

    virtual void PrintMarginalizationSummary() const;

    void UpdateParameterTree();

    unsigned  UpdateFrequency(unsigned N) const;

    std::string fName;

    std::string fSafeName;

    BCParameterSet fParameters;

    BCObservableSet fObservables;

    unsigned fMCMCNChains;

    unsigned fMCMCNLag;

    int fMCMCCurrentIteration;

    unsigned fMCMCNIterationsPreRunCheck;

    unsigned fMCMCPreRunCheckClear;

    int fMCMCNIterationsConvergenceGlobal;

    unsigned fMCMCNIterationsPreRunMax;

    unsigned fMCMCNIterationsRun;

    unsigned fMCMCNIterationsPreRunMin;

    bool fMCMCFlagWriteChainToFile;

    bool fMCMCFlagWritePreRunToFile;

    TFile* fMCMCOutputFile;

    std::string fMCMCOutputFilename;

    std::string fMCMCOutputFileOption;

    double fMCMCScaleFactorLowerLimit;

    double fMCMCScaleFactorUpperLimit;

    std::vector<std::vector<double> > fMCMCProposalFunctionScaleFactor;

    std::vector<double> fMCMCInitialScaleFactors;

    std::vector<TMatrixDSym> fMultivariateProposalFunctionCovariance;

    std::vector<TMatrixD> fMultivariateProposalFunctionCholeskyDecomposition;

    unsigned fMultivariateCovarianceUpdates;

    double fMultivariateCovarianceUpdateLambda;

    double fMultivariateEpsilon;

    double fMultivariateScaleMultiplier;

    std::vector<std::vector<double> > fMCMCInitialPosition;

    double fMCMCEfficiencyMin;

    double fMCMCEfficiencyMax;

    BCEngineMCMC::InitialPositionScheme fInitialPositionScheme;

    unsigned fInitialPositionAttemptLimit;

    bool fMCMCProposeMultivariate;

    double fMCMCProposalFunctionDof;

    BCEngineMCMC::Phase fMCMCPhase;

    std::vector<ChainState> fMCMCStates;

    std::vector<BCEngineMCMC::Statistics> fMCMCStatistics;

    BCEngineMCMC::Statistics fMCMCStatistics_AllChains;

    bool fCorrectRValueForSamplingVariability;

    double fMCMCRValueParametersCriterion;

    std::vector<double> fMCMCRValueParameters;

    TRandom3 fRandom;

    std::vector<TH1*> fH1Marginalized;

    std::vector<std::vector<TH2*> > fH2Marginalized;

    std::vector<std::pair<int, int> > fRequestedH2;

    TTree* fMCMCTree;

    bool fMCMCTreeLoaded;

    bool fMCMCTreeReuseObservables;

    ChainState fMCMCTree_State;

    unsigned int fMCMCTree_Chain;

    TTree* fParameterTree;

    std::vector<double> fLocalModes;

    BCH1D fBCH1DdrawingOptions;

    BCH2D fBCH2DdrawingOptions;

    bool fRescaleHistogramRangesAfterPreRun;

    double fHistogramRescalePadding;

    mutable BCAux::BCTrash<TObject> fObjectTrash;

};

// ---------------------------------------------------------

#endif