BCEngineMCMC.h

#ifndef __BCENGINEMCMC__H
#define __BCENGINEMCMC__H

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

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

#include "BCParameterSet.h"

#include <vector>

class BCH1D;
class BCH2D;

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

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

class BCEngineMCMC
{

   public:

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

      BCEngineMCMC();

      BCEngineMCMC(const BCEngineMCMC & enginemcmc);

      virtual ~BCEngineMCMC();

      BCEngineMCMC & operator = (const BCEngineMCMC & engineMCMC);

      unsigned MCMCGetNChains() const
         { return fMCMCNChains; }

      unsigned MCMCGetNLag() const
         { return fMCMCNLag; }

      const std::vector<unsigned> & MCMCGetNIterations() const
         { return fMCMCNIterations; }

      unsigned MCMCGetCurrentIteration() const
         { return fMCMCCurrentIteration; }

      unsigned MCMCGetCurrentChain() const
         { return fMCMCCurrentChain; }

      unsigned MCMCGetNIterationsConvergenceGlobal() const
         { return fMCMCNIterationsConvergenceGlobal; }

      bool MCMCGetFlagConvergenceGlobal() const
         { return fMCMCFlagConvergenceGlobal; }

      unsigned MCMCGetNIterationsMax() const
         { return fMCMCNIterationsMax; }

      unsigned MCMCGetNIterationsRun() const
         { return fMCMCNIterationsRun; }

      unsigned MCMCGetNIterationsPreRunMin() const
         { return fMCMCNIterationsPreRunMin; }

      unsigned MCMCGetNIterationsUpdate() const
         { return fMCMCNIterationsUpdate; }

      unsigned MCMCGetNIterationsUpdateMax() const
         { return fMCMCNIterationsUpdateMax; }

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

      int MCMCGetNTrials() const
         { return fMCMCNTrials; }

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

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

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

      std::vector<double> MCMCGetTrialFunctionScaleFactor(unsigned ichain) const;

      double MCMCGetTrialFunctionScaleFactor(unsigned ichain, unsigned ipar);

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

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

      double MCMCGetx(unsigned ichain, unsigned ipar) const;

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

      double MCMCGetLogProbx(unsigned ichain);

      int MCMCGetPhase() const
         { return fMCMCPhase; }

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

      std::vector<double> MCMCGetMaximumPoint(unsigned i) const;

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

      int MCMCGetFlagInitialPosition() const
         { return fMCMCFlagInitialPosition; }

      double MCMCGetRValueCriterion() const
         { return fMCMCRValueCriterion; }

      double MCMCGetRValueParametersCriterion() const
         { return fMCMCRValueParametersCriterion; }

      double MCMCGetRValue() const
         { return fMCMCRValue; }

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

      bool MCMCGetRValueStrict() const
         { return fMCMCRValueUseStrict; }

      bool MCMCGetFlagRun() const
         { return fMCMCFlagRun; }

      TTree * MCMCGetMarkovChainTree(unsigned i)
      { if (i < fMCMCTrees.size())
          return fMCMCTrees.at(i);
        else
          return 0; }

      BCH1D * MCMCGetH1Marginalized(unsigned i);

      BCH2D * MCMCGetH2Marginalized(unsigned i, unsigned j);

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

      BCParameter * GetParameter(int index) const
         { return fParameters.Get(index); }

      BCParameter * GetParameter(const char * name) const
         { return fParameters.Get(name); }

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

      unsigned int GetNFixedParameters();

      unsigned int GetNFreeParameters();

      const std::vector<double> & GetBestFitParametersMarginalized() const;

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

      void MCMCSetNChains(unsigned n);

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

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

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

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

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

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

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

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

      void MCMCSetRandomSeed(unsigned seed);

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

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

      void MCMCSetInitialPositions(const 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 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(const std::vector<TTree *> & trees);

      void MCMCInitializeMarkovChainTrees();

      int SetMarginalized(unsigned index, TH1D * h);

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

      void MCMCSetValuesDefault();

      void MCMCSetValuesQuick();

      void MCMCSetValuesDetail();

      void MCMCSetPrecision(BCEngineMCMC::Precision precision);
#if 0

      void SetParameterRange(unsigned int index, double parmin, double parmax);
#endif

      void SetNbins(unsigned int nbins);

      void WriteMarkovChain(bool flag)
      {
         fMCMCFlagWriteChainToFile = flag;
         fMCMCFlagWritePreRunToFile = flag;
      }

#if 0

      void SetMarkovChainInitialPosition(const std::vector<double> & position)
      { fXmetro0 = position; }

      void SetMarkovChainStepSize(double stepsize)
      { fMarkovChainStepSize = stepsize; }

      void SetMarkovChainNIterations(int niterations)
      { fMarkovChainNIterations = niterations;
      fMarkovChainAutoN = false; }

      void SetMarkovChainAutoN(bool flag)
      { fMarkovChainAutoN = flag; }
#endif

      void Copy(const BCEngineMCMC & enginemcmc);

      virtual int AddParameter(const char * name, double min, double max, const char * latexname = "");

      virtual int AddParameter(BCParameter* parameter);

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

      virtual double MCMCTrialFunctionSingle(unsigned ichain, unsigned ipar);

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

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

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

      void MCMCInChainCheckMaximum();

      void MCMCInChainUpdateStatistics();

      void MCMCInChainFillHistograms();

      void MCMCInChainTestConvergenceAllChains();

      void MCMCInChainWriteChains();

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

      int MCMCMetropolis();

      int MCMCMetropolisPreRun();

      void MCMCResetRunStatistics();

      void MCMCInitializeMarkovChains();

      int MCMCInitialize();

      virtual void ResetResults();

      void ClearParameters(bool hard=false)
      { fParameters.Clear(hard); }

      virtual void MCMCIterationInterface();

      virtual void MCMCUserIterationInterface()
         {}

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

   private:


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

      MCMCPointerToGetProposalPoint fMCMCPointerToGetProposalPoint;

      struct MCMCThreadLocalStorage
      {
          std::vector<double> xLocal;

          TRandom3 * rng;

          MCMCThreadLocalStorage(const unsigned & dim);

          MCMCThreadLocalStorage(const MCMCThreadLocalStorage & other);

          MCMCThreadLocalStorage & operator = (const MCMCThreadLocalStorage & other);

          virtual ~MCMCThreadLocalStorage();
      };

      std::vector<MCMCThreadLocalStorage> fMCMCThreadLocalStorage;

      void SyncThreadStorage();

   protected:
      BCParameterSet fParameters;

      unsigned fMCMCNChains;

      unsigned fMCMCNLag;

      std::vector<unsigned> fMCMCNIterations;

      int fMCMCCurrentIteration;

      int fMCMCCurrentChain;

      unsigned fMCMCNIterationsUpdate;

      unsigned fMCMCNIterationsUpdateMax;

      int fMCMCNIterationsConvergenceGlobal;

      bool fMCMCFlagConvergenceGlobal;

      unsigned fMCMCNIterationsMax;

      unsigned fMCMCNIterationsRun;

      unsigned fMCMCNIterationsPreRunMin;

      std::vector<int> fMCMCNTrialsTrue;

      unsigned fMCMCNTrials;

      bool fMCMCFlagWriteChainToFile;

      bool fMCMCFlagWritePreRunToFile;

      std::vector<double> fMCMCTrialFunctionScaleFactor;

      std::vector<double> fMCMCTrialFunctionScaleFactorStart;

      bool fMCMCFlagPreRun;

      bool fMCMCFlagRun;

      std::vector<double> fMCMCInitialPosition;

      std::vector<double> fMCMCEfficiencies;

      double fMCMCEfficiencyMin;

      double fMCMCEfficiencyMax;

      int fMCMCFlagInitialPosition;

      bool fMCMCFlagOrderParameters;

      int fMCMCPhase;

      std::vector<double> fMCMCx;

      std::vector<double> fMCMCxMax;

      std::vector<double> fMCMCxMean;

      std::vector<double> fMCMCxVar;

      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<TH1D *> fMCMCH1Marginalized;
      std::vector<TH2D *> fMCMCH2Marginalized;

      std::vector<TTree *> fMCMCTrees;

      std::vector<double> fMCMCBestFitParameters;

      double fMCMCLogMaximum;

      std::vector<double> fMarginalModes;
};

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

#endif