BCEngineMCMC.h

Go to the documentation of this file.

#ifndef __BCENGINEMCMC__H
#define __BCENGINEMCMC__H

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

#include <vector>

// ROOT classes
class TH1D;
class TH2D;
class TTree;
class TRandom3;

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

class BCEngineMCMC
{

   public:

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

      BCEngineMCMC();

      BCEngineMCMC(int n);
   
      BCEngineMCMC(const BCEngineMCMC & enginemcmc);

      virtual ~BCEngineMCMC();

      BCEngineMCMC & operator = (const BCEngineMCMC & engineMCMC);

      int MCMCGetNParameters()
         { return fMCMCNParameters; }

      int MCMCGetNChains()
         { return fMCMCNChains; }

      int MCMCGetNLag()
         { return fMCMCNLag; }

      std::vector <int> MCMCGetNIterations()
         { return fMCMCNIterations; }

      int MCMCGetCurrentIteration()
         { return fMCMCCurrentIteration; }

      int MCMCGetCurrentChain()
         { return fMCMCCurrentChain; }

      int MCMCGetNIterationsConvergenceGlobal()
         { return fMCMCNIterationsConvergenceGlobal; }

      bool MCMCGetFlagConvergenceGlobal()
         { return fMCMCFlagConvergenceGlobal; }

      int MCMCGetNIterationsMax()
         { return fMCMCNIterationsMax; }

      int MCMCGetNIterationsRun()
         { return fMCMCNIterationsRun; }

      int MCMCGetNIterationsPreRunMin()
         { return fMCMCNIterationsPreRunMin; }

      int MCMCGetNIterationsUpdate()
         { return fMCMCNIterationsUpdate; }

      int MCMCGetNIterationsUpdateMax()
         { return fMCMCNIterationsUpdateMax; }

      std::vector <int> MCMCGetNTrialsTrue()
         { return fMCMCNTrialsTrue; }

      std::vector <int> MCMCGetNTrialsFalse()
         { return fMCMCNTrialsFalse; }

      std::vector <double> MCMCGetprobMean()
         { return fMCMCprobMean; }

      std::vector <double> MCMCGetVariance()
         { return fMCMCprobVar; }

      std::vector <double> MCMCGetTrialFunctionScaleFactor()
         { return fMCMCTrialFunctionScaleFactor; }

      std::vector <double> MCMCGetTrialFunctionScaleFactor(int ichain);

      double MCMCGetTrialFunctionScaleFactor(int ichain, int ipar);

      std::vector <double> MCMCGetx()
         { return fMCMCx; }

      std::vector <double> MCMCGetx(int ichain);

      double MCMCGetx(int ichain, int ipar);

      std::vector <double> MCMCGetLogProbx()
         { return fMCMCprob; }

      double MCMCGetLogProbx(int ichain);

      int MCMCGetPhase()
         { return fMCMCPhase; }

      int MCMCGetCycle()
         { return fMCMCCycle; }

      std::vector <double> MCMCGetMaximumPoints()
         { return fMCMCxMax; }

      std::vector <double> MCMCGetMaximumPoint(int i);

      std::vector <double> MCMCGetMaximumLogProb()
         { return fMCMCprobMax; }

      int MCMCGetFlagInitialPosition()
         { return fMCMCFlagInitialPosition; }

      double MCMCGetRValueCriterion()
         { return fMCMCRValueCriterion; }

      double MCMCGetRValueParametersCriterion()
         { return fMCMCRValueParametersCriterion; }

      double MCMCGetRValue()
         { return fMCMCRValue; }

      double MCMCGetRValueParameters(int i)
         { return fMCMCRValueParameters.at(i); }

      bool MCMCGetRValueStrict()
         { return fMCMCRValueUseStrict; }

      bool MCMCGetFlagRun()
      { return fMCMCFlagRun; }

      TTree * MCMCGetMarkovChainTree(int i)
         { return fMCMCTrees.at(i); }

      TH1D * MCMCGetH1Marginalized(int i);

      TH2D * MCMCGetH2Marginalized(int i, int j);

      TRandom3 * MCMCGetTRandom3()
         { return fRandom; }

      void MCMCSetTrialFunctionScaleFactor(std::vector <double> scale)
         { fMCMCTrialFunctionScaleFactorStart = scale; }

      void MCMCSetNChains(int n);

      void MCMCSetNLag(int n)
         { fMCMCNLag = n; }

      void MCMCSetNIterationsMax(int n)
         { fMCMCNIterationsMax = n; }

      void MCMCSetNIterationsRun(int n)
         { fMCMCNIterationsRun = n; }

      void MCMCSetNIterationsPreRunMin(int n)
         { fMCMCNIterationsPreRunMin = n; }

      void MCMCSetNIterationsUpdate(int n)
         { fMCMCNIterationsUpdate = n; }

      void MCMCSetNIterationsUpdateMax(int n)
         { fMCMCNIterationsUpdateMax = n; }

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

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

      void MCMCSetWriteChainToFile(bool flag)
         { fMCMCFlagWriteChainToFile = flag; }

      void MCMCSetWritePreRunToFile(bool flag)
         { fMCMCFlagWritePreRunToFile = flag; }

      void MCMCSetInitialPositions(std::vector<double> x0s);

      void MCMCSetInitialPositions(std::vector< std::vector<double> > x0s);

      void MCMCSetFlagInitialPosition(int flag)
         { fMCMCFlagInitialPosition = flag; }

      void MCMCSetFlagOrderParameters(bool flag)
         { fMCMCFlagOrderParameters = flag; }

      void MCMCSetFlagFillHistograms(bool flag);

      void MCMCSetFlagFillHistograms(int index, bool flag);

      void MCMCSetFlagPreRun(bool flag) 
      { fMCMCFlagPreRun = flag; }

      void MCMCSetRValueCriterion(double r)
         { fMCMCRValueCriterion = r; }

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

      void MCMCSetRValueStrict(bool strict=true)
      { fMCMCRValueUseStrict = strict; }

      void MCMCSetMarkovChainTrees(std::vector <TTree *> trees);

      void MCMCInitializeMarkovChainTrees();

      int SetMarginalized(int index, TH1D * h);

      int SetMarginalized(int index1, int index2, TH2D * h);

      void MCMCSetValuesDefault();

      void MCMCSetValuesQuick();

      void MCMCSetValuesDetail();

      void MCMCSetPrecision(BCEngineMCMC::Precision precision);


      int MCMCAddParameter(double min, double max);

      virtual void MCMCTrialFunction(int ichain, std::vector <double> &x);

      virtual double MCMCTrialFunctionSingle(int ichain, int ipar);

      bool MCMCGetProposalPointMetropolis(int chain, std::vector <double> &x);

      bool MCMCGetProposalPointMetropolis(int chain, int parameter, std::vector <double> &x);

      bool MCMCGetNewPointMetropolis(int chain = 0);
      bool MCMCGetNewPointMetropolis(int chain, int parameter);

      void MCMCInChainCheckMaximum();

      void MCMCInChainUpdateStatistics();

      void MCMCInChainFillHistograms();

      void MCMCInChainTestConvergenceAllChains();

      double MCMCCalculateRValue(const std::vector<double> & chainMeans,
         const std::vector<double> & chainVariances, double mean, double B,
         double W, int nPoints, int iParam = -1);

      void MCMCInChainWriteChains();

      virtual double LogEval(std::vector <double> parameters);

      int MCMCMetropolis();

      int MCMCMetropolisPreRun();

      void MCMCResetRunStatistics();

      void MCMCInitializeMarkovChains();

      int MCMCInitialize();

      int MCMCResetResults();

      virtual void MCMCIterationInterface()
         {}

      virtual void MCMCCurrentPointInterface(std::vector <double> & point, int ichain, bool accepted)
         {
         //suppress warnings for unused parameters
         //with optimization, no code should be generated
         (void)point;
         (void)ichain;
         (void)accepted;
         }

   private:

      typedef bool (BCEngineMCMC::*MCMCPointerToGetProposalPoint) (int chain, std::vector <double> xnew, std::vector <double> xold) const;

      MCMCPointerToGetProposalPoint fMCMCPointerToGetProposalPoint;

   protected:

      int fMCMCNParameters;

      std::vector <double> fMCMCBoundaryMin;
      std::vector <double> fMCMCBoundaryMax;

      std::vector <bool> fMCMCFlagsFillHistograms;

      int fMCMCNChains;

      int fMCMCNLag;

      std::vector<int> fMCMCNIterations;

      int fMCMCCurrentIteration;

      int fMCMCCurrentChain;

      int fMCMCNIterationsUpdate;

      int fMCMCNIterationsUpdateMax;

      int fMCMCNIterationsConvergenceGlobal;

      bool fMCMCFlagConvergenceGlobal;

      int fMCMCNIterationsMax;

      int fMCMCNIterationsRun;

      int fMCMCNIterationsPreRunMin;

      std::vector<int> fMCMCNTrialsTrue;

      std::vector<int> fMCMCNTrialsFalse;

      bool fMCMCFlagWriteChainToFile;

      bool fMCMCFlagWritePreRunToFile;

      std::vector <double> fMCMCTrialFunctionScaleFactor;


      std::vector <double> fMCMCTrialFunctionScaleFactorStart;

      bool fMCMCFlagPreRun;

      bool fMCMCFlagRun;

      std::vector <double> fMCMCInitialPosition;

      double fMCMCEfficiencyMin;

      double fMCMCEfficiencyMax;

      int fMCMCFlagInitialPosition;

      bool fMCMCFlagOrderParameters;

      bool fMCMCFlagFillHistograms;

      int fMCMCPhase;

      int fMCMCCycle;

      std::vector <double> fMCMCx;

      std::vector <double> fMCMCxMax;

      std::vector <double> fMCMCxMean;

      std::vector <double> fMCMCxVar;

      std::vector <double> fMCMCxLocal;

      std::vector<double> fMCMCprob;

      std::vector <double> fMCMCprobMax;

      std::vector <double> fMCMCprobMean;

      std::vector <double> fMCMCprobVar;

      bool fMCMCRValueUseStrict;

      double fMCMCRValueCriterion;

      double fMCMCRValueParametersCriterion;

      double fMCMCRValue;

      std::vector <double> fMCMCRValueParameters;

      TRandom3 * fRandom;

      std::vector<int> fMCMCH1NBins;

      std::vector <TH1D *> fMCMCH1Marginalized;
      std::vector <TH2D *> fMCMCH2Marginalized;

      std::vector<TTree *> fMCMCTrees;
};

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

#endif