BCIntegrate.h

Go to the documentation of this file.

#ifndef __BCINTEGRATE__H
#define __BCINTEGRATE__H

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

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

#include "BCEngineMCMC.h"

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

namespace BCCubaOptions
{
   struct General
   {
      int ncomp, flags, nregions, neval, fail;
      double error, prob;
      General();
   protected:
      ~General();
   };

   struct Vegas : public General
   {
      int nstart, nincrease, nbatch, gridno;

      Vegas();
   };

   struct Suave : public General
   {
      int nnew;
      double flatness;

      Suave();
   };

   struct Divonne : public General
   {
      int key1, key2, key3, maxpass;
      double border, maxchisq, mindeviation;

      Divonne();
   };

   struct Cuhre : public General
   {
      int key;

      Cuhre();
   };
};

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

class BCIntegrate : public BCEngineMCMC
{

public:

  enum BCOptimizationMethod {
    kOptEmpty,
    kOptSimAnn,
    kOptMetropolis,
    kOptMinuit,
    kOptDefault,
    NOptMethods };

   enum BCIntegrationMethod {
      kIntEmpty,
      kIntMonteCarlo,
      kIntCuba,
      kIntGrid,
      kIntDefault,
      NIntMethods };

   enum BCMarginalizationMethod {
     kMargEmpty,
     kMargMetropolis,
     kMargMonteCarlo,
     kMargGrid,
     kMargDefault,
     NMargMethods };

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

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

   typedef void (BCIntegrate::*tRandomizer)(std::vector<double> &) const;

   typedef double (BCIntegrate::*tEvaluator)(std::vector<double> &, const std::vector<double> &, bool &);

   typedef void (*tIntegralUpdater)(const std::vector<double> &, const int &, double &, double &);

   BCIntegrate();

   BCIntegrate(const BCIntegrate & bcintegrate);

   virtual ~BCIntegrate();

   BCIntegrate & operator = (const BCIntegrate & bcintegrate);

   int ReadMarginalizedFromFile(const char * file);

   BCH1D * GetMarginalized(const BCParameter * parameter);

   BCH1D * GetMarginalized(const char * name)
   { return GetMarginalized(fParameters.Index(name)); }

   BCH1D * GetMarginalized(unsigned index);

   BCH2D * GetMarginalized(const BCParameter * parameter1, const BCParameter * parameter2);

   BCH2D * GetMarginalized(const char * name1, const char * name2)
   { return GetMarginalized(fParameters.Index(name1), fParameters.Index(name2)); }

   BCH2D * GetMarginalized(unsigned index1, unsigned index2);

   int PrintAllMarginalized1D(const char * filebase);
   int PrintAllMarginalized2D(const char * filebase);
   int PrintAllMarginalized(const char * file, std::string options1d="BTsiB3CS1D0pdf0Lmeanmode", std::string options2d="BTfB3CS1meangmode", unsigned int hdiv=1, unsigned int ndiv=1);


   double GetIntegral() const
   { return fIntegral; }

   BCIntegrate::BCOptimizationMethod GetOptimizationMethod() const
     { return fOptimizationMethodCurrent; }

   BCIntegrate::BCIntegrationMethod GetIntegrationMethod() const
   { return fIntegrationMethodCurrent; }

   BCIntegrate::BCMarginalizationMethod GetMarginalizationMethod() const
   { return fMarginalizationMethodCurrent; }

   BCIntegrate::BCSASchedule GetSASchedule() const
   { return fSASchedule; }

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

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

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

   int GetNIterationsMin() const
   { return fNIterationsMin; }

   int GetNIterationsMax() const
   { return fNIterationsMax; }

   int GetNIterationsPrecisionCheck() const
   { return fNIterationsPrecisionCheck; }

   int GetNIterationsOutput() const
   { return fNIterationsOutput; }

   int GetNIterations() const
   { return fNIterations; }

   double GetRelativePrecision() const
   { return fRelativePrecision; }

   double GetAbsolutePrecision() const
   { return fAbsolutePrecision; }

   BCCubaMethod GetCubaIntegrationMethod() const
   { return fCubaIntegrationMethod; }

   const BCCubaOptions::Vegas & GetCubaVegasOptions() const
   { return fCubaVegasOptions; }

   const BCCubaOptions::Suave & GetCubaSuaveOptions() const
   { return fCubaSuaveOptions; }

   const BCCubaOptions::Divonne & GetCubaDivonneOptions() const
   { return fCubaDivonneOptions; }

   const BCCubaOptions::Cuhre & GetCubaCuhreOptions() const
   { return fCubaCuhreOptions; }

   BCH1D* GetSlice(const BCParameter* parameter, const std::vector<double> parameters = std::vector<double>(0), int bins=0, bool flag_norm=true);

   BCH1D* GetSlice(const char * name, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true)
   { return GetSlice(GetParameter(name), parameters, nbins, flag_norm); }

   BCH2D* GetSlice(const BCParameter* parameter1, const BCParameter* parameter2, const std::vector<double> parameters = std::vector<double>(0), int bins=0, bool flag_norm=true);

   BCH2D* GetSlice(const char* name1, const char* name2, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true);

   BCH2D* GetSlice(unsigned index1, unsigned index2, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true);

   double GetError() const
   { return fError; }

   TMinuit * GetMinuit();

   int GetMinuitErrorFlag() const
   { return fMinuitErrorFlag; }

   double GetSAT0() const
   { return fSAT0; }

   double GetSATmin() const
   { return fSATmin; }

   double GetBestFitParameter(unsigned index) const;

   double GetBestFitParameterError(unsigned index) const;

   double GetLogMaximum()
   { return fLogMaximum; };

   const std::vector<double> & GetBestFitParameters() const
   { return fBestFitParameters; }

   const std::vector<double> & GetBestFitParameterErrors() const
   { return fBestFitParameterErrors; }

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

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

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

   void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method);

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

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

   void SetNIterationsMin(int niterations)
   { fNIterationsMin = niterations; }

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

   void SetNIterationsPrecisionCheck(int niterations)
   { fNIterationsPrecisionCheck = niterations; }

   void SetNIterationsOutput(int niterations)
   { fNIterationsOutput = niterations; }

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

   void SetAbsolutePrecision(double absprecision)
   { fAbsolutePrecision = absprecision; }

   void SetCubaIntegrationMethod(BCCubaMethod type);

   void SetCubaOptions(const BCCubaOptions::Vegas & options)
   {  fCubaVegasOptions = options; }

   void SetCubaOptions(const BCCubaOptions::Suave & options)
   {  fCubaSuaveOptions = options; }

   void SetCubaOptions(const BCCubaOptions::Divonne & options)
   {  fCubaDivonneOptions = options; }

   void SetCubaOptions(const BCCubaOptions::Cuhre & options)
   {  fCubaCuhreOptions = options; }

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

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

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

   TTree * GetSATree()
   { return fTreeSA; }

   void InitializeSATree();

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

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

   double Normalize()
   { return Integrate(); };

   double Integrate(BCIntegrationMethod intmethod);

   double Integrate();

   double Integrate(BCIntegrationMethod type, tRandomizer randomizer, tEvaluator evaluator, tIntegralUpdater updater, std::vector<double> &sums);

   // todo document
   double EvaluatorMC(std::vector<double> &sums, const std::vector<double> &point, bool &accepted);
   static void IntegralUpdaterMC(const std::vector<double> &sums, const int &nIterations, double &integral, double &absprecision);

   static int CubaIntegrand(const int * ndim, const double xx[], const int * ncomp, double ff[], void *userdata);
#if 0
   TH1D * Marginalize(BCIntegrationMethod type, unsigned index);

   TH2D * Marginalize(BCIntegrationMethod type, unsigned index1, unsigned index2);

   bool Marginalize(TH1* hist, BCIntegrationMethod type, const std::vector<unsigned> &index);
#endif

   int MarginalizeAll();

   int MarginalizeAll(BCMarginalizationMethod margmethod);

   virtual void MarginalizePreprocess()
   {};

   virtual void MarginalizePostprocess()
   {};

   void SAInitialize();

    std::vector<double> FindMode(std::vector<double> start = std::vector<double>());

    std::vector<double> FindMode(BCIntegrate::BCOptimizationMethod optmethod, std::vector<double> start = std::vector<double>());

   double SATemperature(double t);

   double SATemperatureBoltzmann(double t);

   double SATemperatureCauchy(double t);

   virtual double SATemperatureCustom(double t);

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

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

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

   virtual std::vector<double> GetProposalPointSACustom(const 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 ResetResults();

   std::string DumpIntegrationMethod(BCIntegrationMethod type);

   std::string DumpCurrentIntegrationMethod()
   { return DumpIntegrationMethod(fIntegrationMethodCurrent); }

   std::string DumpUsedIntegrationMethod()
   { return DumpIntegrationMethod(fIntegrationMethodUsed); }

   std::string DumpMarginalizationMethod(BCMarginalizationMethod type);

   std::string DumpCurrentMarginalizationMethod()
   { return DumpMarginalizationMethod(fMarginalizationMethodCurrent); }

   std::string DumpUsedMarginalizationMethod()
   { return DumpMarginalizationMethod(fMarginalizationMethodUsed); }

   std::string DumpOptimizationMethod(BCOptimizationMethod type);

   std::string DumpCurrentOptimizationMethod()
   { return DumpOptimizationMethod(fOptimizationMethodCurrent); }

   std::string DumpUsedOptimizationMethod()
   { return DumpOptimizationMethod(fOptimizationMethodUsed); }

   std::string DumpCubaIntegrationMethod(BCCubaMethod type);

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

   void SetBestFitParameters(const std::vector<double> &x)
   { fBestFitParameters = x; }

   void SetBestFitParameters(const std::vector<double> &x, const double &new_value, double &old_value);

   unsigned GetNIntegrationVariables();

   double CalculateIntegrationVolume();

   bool CheckMarginalizationAvailability(BCMarginalizationMethod type);

   bool CheckMarginalizationIndices(TH1* hist, const std::vector<unsigned> &index);

protected:

   int fID;

   TMinuit * fMinuit;

   double fMinuitArglist[2];
   int fMinuitErrorFlag;

   bool fFlagIgnorePrevOptimization;

   double fSAT0;

   double fSATmin;

   TTree * fTreeSA;

   bool fFlagWriteSAToFile;

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

   std::vector<BCH1D*> fMarginalized1D;

   std::vector<BCH2D*> fMarginalized2D;

 protected:
   unsigned IntegrationOutputFrequency() const;

   void LogOutputAtStartOfIntegration(BCIntegrationMethod type, BCCubaMethod cubatype);
   void LogOutputAtIntegrationStatusUpdate(BCIntegrationMethod type, double integral, double absprecision, int nIterations);
   void LogOutputAtEndOfIntegration(double integral, double absprecision, double relprecision, int nIterations);

   void Copy(const BCIntegrate & bcintegrate);

   bool fFlagMarginalized;

private:

    std::vector<double> FindModeMinuit(std::vector<double> &mode, std::vector<double> &errors, std::vector<double> start = std::vector<double>(0), int printlevel = 1);

   std::vector<double> FindModeMCMC(std::vector<double> &mode, std::vector<double> &errors);

   std::vector<double> FindModeSA(std::vector<double> &mode, std::vector<double> &errors, std::vector<double> start = std::vector<double>(0));

   double IntegrateCuba()
   { return IntegrateCuba(fCubaIntegrationMethod); }

   double IntegrateCuba(BCCubaMethod cubatype);

   double IntegrateSlice();

   BCIntegrate::BCOptimizationMethod fOptimizationMethodCurrent;

   BCIntegrate::BCOptimizationMethod fOptimizationMethodUsed;

   BCIntegrate::BCIntegrationMethod fIntegrationMethodCurrent;

   BCIntegrate::BCIntegrationMethod fIntegrationMethodUsed;

   BCIntegrate::BCMarginalizationMethod fMarginalizationMethodCurrent;

   BCIntegrate::BCMarginalizationMethod fMarginalizationMethodUsed;

   BCIntegrate::BCSASchedule fSASchedule;

   unsigned fNIterationsMin;

   unsigned fNIterationsMax;

   unsigned fNIterationsPrecisionCheck;

   unsigned fNIterationsOutput;

   int fNIterations;

   std::vector<double> fBestFitParameters;

   std::vector<double> fBestFitParameterErrors;

   double fLogMaximum;

   double fIntegral;

   double fRelativePrecision;

   double fAbsolutePrecision;

   double fError;

   BCCubaMethod fCubaIntegrationMethod;
   BCCubaOptions::Vegas fCubaVegasOptions;
   BCCubaOptions::Suave fCubaSuaveOptions;
   BCCubaOptions::Divonne fCubaDivonneOptions;
   BCCubaOptions::Cuhre fCubaCuhreOptions;

};

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

#endif