BCIntegrate.h

Go to the documentation of this file.

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

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

#ifndef __BCINTEGRATE__H
#define __BCINTEGRATE__H

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

#include <iostream>
#include <vector>

#include <math.h>

#include <TROOT.h>
#include <TRandom3.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TMinuit.h>
#include <TTree.h>

#include "BCMath.h"
#include "BCParameter.h"
#include "BCDataPoint.h"
#include "BCEngineMCMC.h"

#define DEBUG 0

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

class BCIntegrate : public BCEngineMCMC
{

 public:   

   /* @{ */ 

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

   enum BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis }; 

   enum BCOptimizationMethod { kOptSimulatedAnnealing, kOptMetropolis, kOptMinuit }; 

   /* @} */ 

   /* @{ */ 

   BCIntegrate();

   BCIntegrate(BCParameterSet * par);

   virtual ~BCIntegrate();

   /* @} */ 

   /* @{ */ 

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

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

   BCIntegrate::BCOptimizationMethod GetOptimizationMethod()
      { return fOptimizationMethod; };

   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 GetError()
   { return fError; }; 

   int GetNbins()
   { return fNbins; };

   TMinuit * GetMinuit(); 

   int GetMinuitErrorFlag() 
   { return fMinuitErrorFlag; }; 

   TTree * GetMarkovChainTree()
      { return fMarkovChainTree; }; 

   void GetRandomPointSA(std::vector <double> &x, double T, double step);

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

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

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

   /* @} */ 

   /* @{ */ 

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

   void SetParameters(BCParameterSet * par);

   void SetVarList(int * varlist);

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

   void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method)
   { fIntegrationMethod = method; };

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

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

   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 n);
   void SetNbinsX(int n);
   void SetNbinsY(int n);

   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 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.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 FindMode();

   void FindModeSA();

   void FindModeMinuit(); 

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

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

   /* @} */ 

 private:

   BCParameterSet * fx;

   double * fMin;

   double * fMax;

   int * fVarlist;

   int fNiterPerDimension;

   BCIntegrate::BCIntegrationMethod fIntegrationMethod;

   BCIntegrate::BCMarginalizationMethod fMarginalizationMethod;
  
   BCIntegrate::BCOptimizationMethod fOptimizationMethod;

   int fNIterationsMax; 

   int fNIterations; 

   double fRelativePrecision;

   double fError; 

   int fNmetro;
   int fNacceptedMCMC; 

   std::vector <double> fXmetro0;

   std::vector <double> fXmetro1;

   double fMarkovChainValue; 

   /*
    * Interface to BCEngineMCMC::MCMCUpdateFunctionFitting.
    */
   void MCMCUpdateFunctionFitting();

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

   std::vector <TH1D *> fHProb1D;

   std::vector <TH2D *> fHProb2D;

   int fFitFunctionIndexX; 
   int fFitFunctionIndexY; 

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

   TH2D * fErrorBandXY; 
   int fErrorBandNbinsX; 
   int fErrorBandNbinsY; 

   TMinuit * fMinuit; 

   double fMinuitArglist[2]; 
   int fMinuitErrorFlag; 

   bool fFlagWriteMarkovChain; 

   TTree * fMarkovChainTree; 

   int fMCMCIteration; 

};

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

#endif

---