BCIntegrate.h

Go to the documentation of this file.

#ifndef __BCINTEGRATE__H
#define __BCINTEGRATE__H

/*
 * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger.
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

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

#include <vector>

#include <math.h>

#include "BAT/BCEngineMCMC.h"
#include "BAT/BCParameter.h"
#include "BAT/BCDataPoint.h"

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


#define BC_DEBUG 0

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

class BCIntegrate : public BCEngineMCMC
{

   public:

      /* @{ */

      enum BCIntegrationMethod { kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba };

      enum BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis };

      enum BCOptimizationMethod { kOptSA, kOptMetropolis, kOptMinuit };
      
      enum BCSASchedule { kSACauchy, kSABoltzmann, kSACustom };

      /* @} */

      /* @{ */

      BCIntegrate();

      BCIntegrate(BCParameterSet * par);

      virtual ~BCIntegrate();

      /* @} */

      /* @{ */

      BCIntegrate::BCIntegrationMethod GetIntegrationMethod()
         { return fIntegrationMethod; };

      BCIntegrate::BCMarginalizationMethod GetMarginalizationMethod()
         { return fMarginalizationMethod; };

      BCIntegrate::BCOptimizationMethod GetOptimizationMethod()
         { return fOptimizationMethod; };
      
      BCIntegrate::BCOptimizationMethod GetOptimizationMethodMode()
         { return fOptimizationMethodMode; };
      
      BCIntegrate::BCSASchedule GetSASchedule()
         { return fSASchedule; };

      void GetRandomVector(std::vector <double> &x);

      virtual void GetRandomVectorMetro(std::vector <double> &x);

      double GetRandomPoint(std::vector <double> &x);

      double GetRandomPointImportance(std::vector <double> &x);

      void GetRandomPointMetro(std::vector <double> &x);

      void GetRandomPointSamplingMetro(std::vector <double> &x);

      int GetNiterationsPerDimension()
         { return fNiterPerDimension; };

      int GetNSamplesPer2DBin()
         { return fNSamplesPer2DBin; };

      int GetNvar()
         { return fNvar; };

      int GetNIterationsMax()
         { return fNIterationsMax; };

      int GetNIterations()
         { return fNIterations; };

      double GetRelativePrecision()
         { return fRelativePrecision; };

      /*
       * @return The uncertainty in the most recent Monte Carlo integration */
      double GetError()
         { return fError; };

      /*
       * @return number of bins per dimension for the marginalized distributions */
      int GetNbins()
         { return fNbins; };

      /*
       * @return Minuit used for mode finding */
      TMinuit * GetMinuit();

      int GetMinuitErrorFlag()
         { return fMinuitErrorFlag; };

      /*
       * @return The ROOT tree containing the Markov chain */
      TTree * GetMarkovChainTree()
         { return fMarkovChainTree; };

      /*
       * Returns the actual point in the markov chain */
      std::vector<double> * GetMarkovChainPoint()
         { return &fXmetro1; };

      int * GetMCMCIteration()
         { return &fMCMCIteration; };

      double * GetMarkovChainValue()
         { return &fMarkovChainValue; };

      double GetSAT0()
         { return fSAT0; }

      double GetSATmin()
         { return fSATmin; }

      double * GetSAP2Temperature()
         { return &fSATemperature; };

      double * GetSAP2LogProb()
         { return &fSALogProb; }; 
      
      std::vector<double> * GetSAP2x()
         { return &fSAx; }; 

      int * GetSAP2NIterations()
         { return &fSANIterations; };

      /* @} */

      /* @{ */

      void SetMinuitArlist(double * arglist)
         { fMinuitArglist[0] = arglist[0];
           fMinuitArglist[1] = arglist[1]; };

      void SetFlagIgnorePrevOptimization(bool flag)
      { fFlagIgnorePrevOptimization = flag; }; 

      void SetParameters(BCParameterSet * par);

      void SetVarList(int * varlist);

      void SetVar(int index){fVarlist[index]=1;};

      void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method);

      void SetMarginalizationMethod(BCIntegrate::BCMarginalizationMethod method)
         { fMarginalizationMethod = method; };

      void SetOptimizationMethod(BCIntegrate::BCOptimizationMethod method)
         { fOptimizationMethod = method; };

      void SetSASchedule(BCIntegrate::BCSASchedule schedule)
         { fSASchedule = schedule; };

      void SetNiterationsPerDimension(int niterations)
         { fNiterPerDimension = niterations; };

      void SetNSamplesPer2DBin(int n)
         { fNSamplesPer2DBin = n; };

      void SetNIterationsMax(int niterations)
         { fNIterationsMax = niterations; };

      void SetRelativePrecision(double relprecision)
         { fRelativePrecision = relprecision; };

      void SetNbins(int nbins, int index = -1);

//    void SetNbins(int n);
//    void SetNbinsX(int n);
//    void SetNbinsY(int n);

      /*
       * Turn on or off the filling of the error band during the MCMC run.
       * @param flag set to true for turning on the filling */
      void SetFillErrorBand(bool flag = true)
         { fFillErrorBand=flag; };

      /*
       * Turn off filling of the error band during the MCMC run.
       * This method is equivalent to SetFillErrorBand(false) */
      void UnsetFillErrorBand()
         { fFillErrorBand=false; };

      void SetFitFunctionIndexX(int index)
         { fFitFunctionIndexX = index; };

      void SetFitFunctionIndexY(int index)
         { fFitFunctionIndexY = index; };

      void SetFitFunctionIndices(int indexx, int indexy)
         { this -> SetFitFunctionIndexX(indexx);
            this -> SetFitFunctionIndexY(indexy); };

      void SetDataPointLowerBoundaries(BCDataPoint * datasetlowerboundaries)
         { fDataPointLowerBoundaries = datasetlowerboundaries; };

      void SetDataPointUpperBoundaries(BCDataPoint * datasetupperboundaries)
         { fDataPointUpperBoundaries = datasetupperboundaries; };

      void SetDataPointLowerBoundary(int index, double lowerboundary)
         { fDataPointLowerBoundaries -> SetValue(index, lowerboundary); };

      void SetDataPointUpperBoundary(int index, double upperboundary)
         { fDataPointUpperBoundaries -> SetValue(index, upperboundary); };

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

      void SetMarkovChainTree(TTree * tree)
         { fMarkovChainTree = tree; };

      void SetMarkovChainInitialPosition(std::vector<double> position);

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

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

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

      void SetMode(std::vector <double> mode);

      void SetErrorBandHisto(TH2D * h)
         { fErrorBandXY = h; };

      void SetSAT0(double T0)
         { fSAT0 = T0; }

      void SetSATmin(double Tmin)
         { fSATmin = Tmin; }

      void SetFlagWriteSAToFile(bool flag)
      { fFlagWriteSAToFile = flag; }; 

      /*
       * Sets the tree containing the Markov chains. */
      void SetSATree(TTree* tree)
      { fTreeSA = tree; };

      /* @} */

      /* @{ */

      void DeleteVarList();

      void ResetVarlist(int v);

      void UnsetVar(int index)
         { fVarlist[index] = 0; };

      virtual double Eval(std::vector <double> x);

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

      virtual double EvalSampling(std::vector <double> x);

      double LogEvalSampling(std::vector <double> x);

      double EvalPrint(std::vector <double> x);

      virtual double FitFunction(std::vector <double> x, std::vector <double> parameters)
         { return 0.; };

      double Integrate();

      double IntegralMC(std::vector <double> x, int * varlist);

      double IntegralMC(std::vector <double> x);

      double IntegralMetro(std::vector <double> x);

      double IntegralImportance(std::vector <double> x);

      double CubaIntegrate(int method, std::vector<double> parameters_double, std::vector<double> parameters_int);

      double CubaIntegrate();

      static void CubaIntegrand(const int * ndim, const double xx[], const int * ncomp, double ff[]);

      TH1D * Marginalize(BCParameter * parameter);

      TH2D * Marginalize(BCParameter * parameter1, BCParameter * parameter2);

      TH1D * MarginalizeByIntegrate(BCParameter * parameter);

      TH2D * MarginalizeByIntegrate(BCParameter * parameter1, BCParameter * parameter2);

      TH1D * MarginalizeByMetro(BCParameter * parameter);

      TH2D * MarginalizeByMetro(BCParameter * parameter1, BCParameter * parameter2);

      int MarginalizeAllByMetro(const char * name);

      TH1D * GetH1D(int parIndex);

      int GetH2DIndex(int parIndex1, int parIndex2);

      TH2D * GetH2D(int parIndex1, int parIndex2);

      void InitMetro();

      void SAInitialize();

//    void FindMode();

      virtual void FindModeMinuit(std::vector<double> start = std::vector<double>(0), int printlevel = 1);

      void FindModeMCMC();
//    void FindModeMCMC(int flag_run = 0);

      void FindModeSA(std::vector<double> start = std::vector<double>(0));

      double SATemperature(double t);

      double SATemperatureBoltzmann(double t);

      double SATemperatureCauchy(double t);
      
      virtual double SATemperatureCustom(double t);

      std::vector<double> GetProposalPointSA(std::vector<double> x, int t);

      std::vector<double> GetProposalPointSABoltzmann(std::vector<double> x, int t);

      std::vector<double> GetProposalPointSACauchy(std::vector<double> x, int t);

      virtual std::vector<double> GetProposalPointSACustom(std::vector<double> x, int t);

      std::vector<double> SAHelperGetRandomPointOnHypersphere();

      double SAHelperGetRadialCauchy();

      double SAHelperSinusToNIntegral(int dim, double theta);
      

      static void FCNLikelihood(int &npar, double * grad, double &fval, double * par, int flag);

      /*
       * Method executed for every iteration of the MCMC. User's code should be
       * provided via overloading in the derived class*/
      virtual void MCMCUserIterationInterface()
         {};

      /* @} */

   private:

      BCParameterSet * fx;

      double * fMin;

      double * fMax;

      int * fVarlist;

      int fNiterPerDimension;

      BCIntegrate::BCIntegrationMethod fIntegrationMethod;

      BCIntegrate::BCMarginalizationMethod fMarginalizationMethod;

      BCIntegrate::BCOptimizationMethod fOptimizationMethod;

      BCIntegrate::BCOptimizationMethod fOptimizationMethodMode; 

      BCIntegrate::BCSASchedule fSASchedule;

      int fNIterationsMax;

      int fNIterations;

      double fRelativePrecision;

      double fError;

      int fNmetro;
      int fNacceptedMCMC;

      std::vector <double> fXmetro0;

      std::vector <double> fXmetro1;

      double fMarkovChainValue;

      /*
       * Method executed for every iteration of the MCMC, overloaded from BCEngineMCMC. */
      void MCMCIterationInterface();

      /*
       * Fill error band histogram for curreent iteration. This method is called from MCMCIterationInterface() */
      void MCMCFillErrorBand();

   protected:

      int fNvar;

      int fNbins;

      int fNSamplesPer2DBin;

      double fMarkovChainStepSize;

      int fMarkovChainNIterations;

      bool fMarkovChainAutoN;

      BCDataPoint * fDataPointLowerBoundaries;

      BCDataPoint * fDataPointUpperBoundaries;

      std::vector <bool> fDataFixedValues;

      TRandom3 * fRandom;

      std::vector <double> fBestFitParameters;
      std::vector <double> fBestFitParameterErrors;

      std::vector <double> fBestFitParametersMarginalized;

      std::vector <TH1D *> fHProb1D;

      std::vector <TH2D *> fHProb2D;

      bool fFillErrorBand;

      int fFitFunctionIndexX;
      int fFitFunctionIndexY;

      bool fErrorBandContinuous;
      std::vector <double> fErrorBandX;

      TH2D * fErrorBandXY;

      int fErrorBandNbinsX;

      int fErrorBandNbinsY;

      TMinuit * fMinuit;

      double fMinuitArglist[2];
      int fMinuitErrorFlag;

      double fFlagIgnorePrevOptimization; 

      bool fFlagWriteMarkovChain;

      TTree * fMarkovChainTree;

      int fMCMCIteration;

      double fSAT0;

      double fSATmin;

      TTree * fTreeSA; 

      bool fFlagWriteSAToFile; 

      int fSANIterations; 
      double fSATemperature; 
      double fSALogProb; 
      std::vector<double> fSAx; 

};

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

#endif

---