BCIntegrate.h

Go to the documentation of this file.

#ifndef __BCINTEGRATE__H
#define __BCINTEGRATE__H

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

#include <vector>

#include <math.h>

#include "BCEngineMCMC.h"
#include "BCParameter.h"
#include "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, NIntMethods };

      enum BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis, NMargMethods };

      enum BCOptimizationMethod { kOptSA, kOptMetropolis, kOptMinuit, NOptMethods };

      enum BCSASchedule { kSACauchy, kSABoltzmann, kSACustom, NSAMethods };

      enum BCCubaMethod { kCubaVegas, kCubaSuave, kCubaDivonne, kCubaCuhre };

      BCIntegrate();

      BCIntegrate(BCParameterSet * par);

      BCIntegrate(const BCIntegrate & bcintegrate);

      virtual ~BCIntegrate();

         BCIntegrate & operator = (const BCIntegrate & 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; }

      double GetAbsolutePrecision()
         { return fAbsolutePrecision; }

      BCCubaMethod GetCubaIntegrationMethod()
         { return fCubaIntegrationMethod; }

      int GetCubaMinEval()
         { return fCubaMinEval; }

      int GetCubaMaxEval()
         { return fCubaMaxEval; }

      int GetCubaVerbositylevel()
         { return fCubaVerbosity; }

      int GetCubaVegasNStart()
         { return fCubaVegasNStart; }

      int GetCubaVegasNIncrease()
         { return fCubaVegasNIncrease; }

      int GetCubaSuaveNNew()
         { return fCubaSuaveNNew; }

      double GetCubaSuaveFlatness()
         { return fCubaSuaveFlatness; }

      double GetError()
         { return fError; }

      int GetNbins()
         { return fNbins; }

      TMinuit * GetMinuit();

      int GetMinuitErrorFlag()
         { return fMinuitErrorFlag; }

      TTree * GetMarkovChainTree()
         { return fMarkovChainTree; }

      std::vector<double> * GetMarkovChainPoint()
         { return &fXmetro1; }

      int * GetMCMCIteration()
         { return &fMCMCIteration; }

      double * GetMarkovChainValue()
         { return &fMarkovChainValue; }

      double GetSAT0()
         { return fSAT0; }

      double GetSATmin()
         { return fSATmin; }

      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 SetOptimizationMethodMode(BCIntegrate::BCOptimizationMethod method)
         { fOptimizationMethodMode = 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 SetAbsolutePrecision(double absprecision)
         { fAbsolutePrecision = absprecision; }

      void SetCubaIntegrationMethod(BCCubaMethod type);

      void SetCubaMinEval(int n)
         { fCubaMinEval = n; }

      void SetCubaMaxEval(int n)
         { fCubaMaxEval = n; }

      void SetCubaVerbosityLevel(int n)
         { fCubaVerbosity = n; }

      void SetCubaVegasNStart(int n)
         { fCubaVegasNStart = n; }

      void SetCubaVegasNIncrease(int n)
         { fCubaVegasNIncrease = n; }

      void SetCubaSuaveNNew(int n)
         { fCubaSuaveNNew = n; }

      void SetCubaSuaveFlatness(double p)
         { fCubaSuaveFlatness = p; }

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

      void SetFillErrorBand(bool flag = true)
         { fFillErrorBand=flag; }

      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;
           fMCMCFlagWritePreRunToFile = 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; }

      void SetSATree(TTree * tree)
         { fTreeSA = tree; }

      TTree * GetSATree()
         { return fTreeSA; }

      void InitializeSATree();

      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(BCIntegrate::BCCubaMethod 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[]);
      static int CubaIntegrand(const int * ndim, const double xx[], const int * ncomp, double ff[], void *userdata);

      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);

      virtual void MCMCUserIterationInterface()
         {}

      int IntegrateResetResults();

      std::string DumpIntegrationMethod(BCIntegrationMethod type);
      std::string DumpIntegrationMethod()
         { return DumpIntegrationMethod(fIntegrationMethod); }

      std::string DumpMarginalizationMethod(BCMarginalizationMethod type);
      std::string DumpMarginalizationMethod()
         { return DumpMarginalizationMethod(fMarginalizationMethod); }

      std::string DumpOptimizationMethod(BCOptimizationMethod type);
      std::string DumpOptimizationMethod()
         { return DumpOptimizationMethod(fOptimizationMethod); }
      std::string DumpUsedOptimizationMethod()
         { return DumpOptimizationMethod(fOptimizationMethodMode); }

      std::string DumpCubaIntegrationMethod(BCCubaMethod type);
      std::string DumpCubaIntegrationMethod()
         { return DumpCubaIntegrationMethod(fCubaIntegrationMethod); }


   protected:

      int fNvar;

      int fNbins;

      int fNSamplesPer2DBin;

      double fMarkovChainStepSize;

      int fMarkovChainNIterations;

      bool fMarkovChainAutoN;

      BCDataPoint * fDataPointLowerBoundaries;

      BCDataPoint * fDataPointUpperBoundaries;

      std::vector <bool> fDataFixedValues;

      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;

   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 fAbsolutePrecision;

      BCCubaMethod fCubaIntegrationMethod;

      int fCubaMinEval;

      int fCubaMaxEval;

      int fCubaVerbosity;

      int fCubaVegasNStart;

      int fCubaVegasNIncrease;

      int fCubaSuaveNNew;

      double fCubaSuaveFlatness;

      double fError;

      int fNmetro;
      int fNacceptedMCMC;

      std::vector <double> fXmetro0;

      std::vector <double> fXmetro1;

      double fMarkovChainValue;

      void MCMCIterationInterface();

      void MCMCFillErrorBand();
};

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

#endif