BAT  0.9.4
The Bayesian analysis toolkit
 All Classes Namespaces Functions Variables Enumerations
BCIntegrate.h
1 #ifndef __BCINTEGRATE__H
2 #define __BCINTEGRATE__H
3 
16 /*
17  * Copyright (C) 2007-2014, the BAT core developer team
18  * All rights reserved.
19  *
20  * For the licensing terms see doc/COPYING.
21  * For documentation see http://mpp.mpg.de/bat
22  */
23 
24 // ---------------------------------------------------------
25 
26 #include "BCEngineMCMC.h"
27 
28 // ROOT classes
29 class TH1;
30 class TH1D;
31 class TH2D;
32 class TMinuit;
33 class TTree;
34 
40 namespace BCCubaOptions
41 {
42 struct General
43 {
44  int ncomp, flags, nregions, neval, fail;
45  double error, prob;
46  General();
47 protected:
48  ~General();
49 };
50 
51 struct Vegas : public General
52 {
53  int nstart, nincrease, nbatch, gridno;
54 
55  Vegas();
56 };
57 
58 struct Suave : public General
59 {
60  int nnew;
61  double flatness;
62 
63  Suave();
64 };
65 
66 struct Divonne : public General
67 {
68  int key1, key2, key3, maxpass;
69  double border, maxchisq, mindeviation;
70 
71  Divonne();
72 };
73 
74 struct Cuhre : public General
75 {
76  int key;
77 
78  Cuhre();
79 };
80 }
81 // ---------------------------------------------------------
82 
83 class BCIntegrate : public BCEngineMCMC
84 {
85 
86 public:
87 
93  enum BCOptimizationMethod {
94  kOptEmpty,
95  kOptSimAnn,
96  kOptMetropolis,
97  kOptMinuit,
98  kOptDefault,
99  NOptMethods };
100 
103  enum BCIntegrationMethod {
104  kIntEmpty,
105  kIntMonteCarlo,
106  kIntCuba,
107  kIntGrid,
108  kIntDefault,
109  NIntMethods };
110 
113  enum BCMarginalizationMethod {
114  kMargEmpty,
115  kMargMetropolis,
116  kMargMonteCarlo,
117  kMargGrid,
118  kMargDefault,
119  NMargMethods };
120 
123  enum BCSASchedule { kSACauchy, kSABoltzmann, kSACustom, NSAMethods };
124 
127  enum BCCubaMethod { kCubaVegas, kCubaSuave, kCubaDivonne, kCubaCuhre, NCubaMethods};
128 
136  typedef void (BCIntegrate::*tRandomizer)(std::vector<double> &) const;
137 
140  typedef double (BCIntegrate::*tEvaluator)(std::vector<double> &, const std::vector<double> &, bool &);
141 
144  typedef void (*tIntegralUpdater)(const std::vector<double> &, const int &, double &, double &);
145 
154  BCIntegrate();
155 
158  BCIntegrate(const BCIntegrate & bcintegrate);
159 
162  virtual ~BCIntegrate();
163 
170  BCIntegrate & operator = (const BCIntegrate & bcintegrate);
171 
178  int ReadMarginalizedFromFile(const char * file);
179 
187  BCH1D * GetMarginalized(const BCParameter * parameter);
188 
196  BCH1D * GetMarginalized(const char * name)
197  { return GetMarginalized(fParameters.Index(name)); }
198 
206  BCH1D * GetMarginalized(unsigned index);
207 
217  BCH2D * GetMarginalized(const BCParameter * parameter1, const BCParameter * parameter2);
218 
227  BCH2D * GetMarginalized(const char * name1, const char * name2)
228  { return GetMarginalized(fParameters.Index(name1), fParameters.Index(name2)); }
229 
238  BCH2D * GetMarginalized(unsigned index1, unsigned index2);
239 
242  int PrintAllMarginalized1D(const char * filebase);
243  int PrintAllMarginalized2D(const char * filebase);
244  int PrintAllMarginalized(const char * file, std::string options1d="BTsiB3CS1D0pdf0Lmeanmode", std::string options2d="BTfB3CS1meangmode", unsigned int hdiv=1, unsigned int ndiv=1);
245 
246 
249  double GetIntegral() const
250  { return fIntegral; }
251 
254  BCIntegrate::BCOptimizationMethod GetOptimizationMethod() const
255  { return fOptimizationMethodCurrent; }
256 
259  BCIntegrate::BCIntegrationMethod GetIntegrationMethod() const
260  { return fIntegrationMethodCurrent; }
261 
264  BCIntegrate::BCMarginalizationMethod GetMarginalizationMethod() const
265  { return fMarginalizationMethodCurrent; }
266 
269  BCIntegrate::BCSASchedule GetSASchedule() const
270  { return fSASchedule; }
271 
275  void GetRandomVectorUnitHypercube(std::vector<double> &x) const;
276 
280  void GetRandomVectorInParameterSpace(std::vector<double> &x) const;
281 
287  double GetRandomPoint(std::vector<double> &x);
288 
291  int GetNIterationsMin() const
292  { return fNIterationsMin; }
293 
296  int GetNIterationsMax() const
297  { return fNIterationsMax; }
298 
301  int GetNIterationsPrecisionCheck() const
302  { return fNIterationsPrecisionCheck; }
303 
306  int GetNIterationsOutput() const
307  { return fNIterationsOutput; }
308 
311  int GetNIterations() const
312  { return fNIterations; }
313 
316  double GetRelativePrecision() const
317  { return fRelativePrecision; }
318 
321  double GetAbsolutePrecision() const
322  { return fAbsolutePrecision; }
323 
326  BCCubaMethod GetCubaIntegrationMethod() const
327  { return fCubaIntegrationMethod; }
328 
332  const BCCubaOptions::Vegas & GetCubaVegasOptions() const
333  { return fCubaVegasOptions; }
334 
335  const BCCubaOptions::Suave & GetCubaSuaveOptions() const
336  { return fCubaSuaveOptions; }
337 
338  const BCCubaOptions::Divonne & GetCubaDivonneOptions() const
339  { return fCubaDivonneOptions; }
340 
341  const BCCubaOptions::Cuhre & GetCubaCuhreOptions() const
342  { return fCubaCuhreOptions; }
343 
351  BCH1D* GetSlice(const BCParameter* parameter, const std::vector<double> parameters = std::vector<double>(0), int bins=0, bool flag_norm=true);
352 
360  BCH1D* GetSlice(const char * name, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true)
361  { return GetSlice(GetParameter(name), parameters, nbins, flag_norm); }
362 
371  BCH2D* GetSlice(const BCParameter* parameter1, const BCParameter* parameter2, const std::vector<double> parameters = std::vector<double>(0), int bins=0, bool flag_norm=true);
372 
381  BCH2D* GetSlice(const char* name1, const char* name2, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true);
382 
391  BCH2D* GetSlice(unsigned index1, unsigned index2, const std::vector<double> parameters = std::vector<double>(0), int nbins=0, bool flag_norm=true);
392 
395  double GetError() const
396  { return fError; }
397 
400  TMinuit * GetMinuit();
401 
404  int GetMinuitErrorFlag() const
405  { return fMinuitErrorFlag; }
406 
409  double GetSAT0() const
410  { return fSAT0; }
411 
414  double GetSATmin() const
415  { return fSATmin; }
416 
422  double GetBestFitParameter(unsigned index) const;
423 
429  double GetBestFitParameterError(unsigned index) const;
430 
434  double GetLogMaximum()
435  { return fLogMaximum; };
436 
441  const std::vector<double> & GetBestFitParameters() const
442  { return fBestFitParameters; }
443 
446  const std::vector<double> & GetBestFitParameterErrors() const
447  { return fBestFitParameterErrors; }
448 
456  void SetMinuitArlist(double * arglist)
457  { fMinuitArglist[0] = arglist[0];
458  fMinuitArglist[1] = arglist[1]; }
459 
462  void SetFlagIgnorePrevOptimization(bool flag)
463  { fFlagIgnorePrevOptimization = flag; }
464 
467  void SetOptimizationMethod(BCIntegrate::BCOptimizationMethod method)
468  { fOptimizationMethodCurrent = method; }
469 
472  void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method);
473 
476  void SetMarginalizationMethod(BCIntegrate::BCMarginalizationMethod method)
477  { fMarginalizationMethodCurrent = method; }
478 
481  void SetSASchedule(BCIntegrate::BCSASchedule schedule)
482  { fSASchedule = schedule; }
483 
486  void SetNIterationsMin(int niterations)
487  { fNIterationsMin = niterations; }
488 
491  void SetNIterationsMax(int niterations)
492  { fNIterationsMax = niterations; }
493 
496  void SetNIterationsPrecisionCheck(int niterations)
497  { fNIterationsPrecisionCheck = niterations; }
498 
501  void SetNIterationsOutput(int niterations)
502  { fNIterationsOutput = niterations; }
503 
507  void SetRelativePrecision(double relprecision)
508  { fRelativePrecision = relprecision; }
509 
512  void SetAbsolutePrecision(double absprecision)
513  { fAbsolutePrecision = absprecision; }
514 
517  void SetCubaIntegrationMethod(BCCubaMethod type);
518 
522  void SetCubaOptions(const BCCubaOptions::Vegas & options)
523  { fCubaVegasOptions = options; }
524 
525  void SetCubaOptions(const BCCubaOptions::Suave & options)
526  { fCubaSuaveOptions = options; }
527 
528  void SetCubaOptions(const BCCubaOptions::Divonne & options)
529  { fCubaDivonneOptions = options; }
530 
531  void SetCubaOptions(const BCCubaOptions::Cuhre & options)
532  { fCubaCuhreOptions = options; }
533 
536  void SetSAT0(double T0)
537  { fSAT0 = T0; }
538 
541  void SetSATmin(double Tmin)
542  { fSATmin = Tmin; }
543 
544  void SetFlagWriteSAToFile(bool flag)
545  { fFlagWriteSAToFile = flag; }
546 
549  TTree * GetSATree()
550  { return fTreeSA; }
551 
554  void InitializeSATree();
555 
566  virtual double Eval(const std::vector<double> &x);
567 
573  virtual double LogEval(const std::vector<double> &x);
574 
577  double Normalize()
578  { return Integrate(); };
579 
584  double Integrate(BCIntegrationMethod intmethod);
585 
590  double Integrate();
591 
601  double Integrate(BCIntegrationMethod type, tRandomizer randomizer, tEvaluator evaluator, tIntegralUpdater updater, std::vector<double> &sums);
602 
603  // todo document
604  double EvaluatorMC(std::vector<double> &sums, const std::vector<double> &point, bool &accepted);
605  static void IntegralUpdaterMC(const std::vector<double> &sums, const int &nIterations, double &integral, double &absprecision);
606 
614  static int CubaIntegrand(const int * ndim, const double xx[], const int * ncomp, double ff[], void *userdata);
615 #if 0
616  TH1D * Marginalize(BCIntegrationMethod type, unsigned index);
617 
618  TH2D * Marginalize(BCIntegrationMethod type, unsigned index1, unsigned index2);
619 
620  bool Marginalize(TH1* hist, BCIntegrationMethod type, const std::vector<unsigned> &index);
621 #endif
622 
627  int MarginalizeAll();
628 
635  int MarginalizeAll(BCMarginalizationMethod margmethod);
636 
640  virtual void MarginalizePreprocess()
641  {};
642 
646  virtual void MarginalizePostprocess()
647  {};
648 
651  void SAInitialize();
652 
664  std::vector<double> FindMode(std::vector<double> start = std::vector<double>());
665 
672  std::vector<double> FindMode(BCIntegrate::BCOptimizationMethod optmethod, std::vector<double> start = std::vector<double>());
673 
679  double SATemperature(double t);
680 
685  double SATemperatureBoltzmann(double t);
686 
691  double SATemperatureCauchy(double t);
692 
698  virtual double SATemperatureCustom(double t);
699 
707  std::vector<double> GetProposalPointSA(const std::vector<double> &x, int t);
708 
715  std::vector<double> GetProposalPointSABoltzmann(const std::vector<double> &x, int t);
716 
723  std::vector<double> GetProposalPointSACauchy(const std::vector<double> &x, int t);
724 
732  virtual std::vector<double> GetProposalPointSACustom(const std::vector<double> &x, int t);
733 
738  std::vector<double> SAHelperGetRandomPointOnHypersphere();
739 
743  double SAHelperGetRadialCauchy();
744 
748  double SAHelperSinusToNIntegral(int dim, double theta);
749 
750 
751  static void FCNLikelihood(int &npar, double * grad, double &fval, double * par, int flag);
752 
755  virtual void ResetResults();
756 
761  std::string DumpIntegrationMethod(BCIntegrationMethod type);
762 
766  std::string DumpCurrentIntegrationMethod()
767  { return DumpIntegrationMethod(fIntegrationMethodCurrent); }
768 
772  std::string DumpUsedIntegrationMethod()
773  { return DumpIntegrationMethod(fIntegrationMethodUsed); }
774 
779  std::string DumpMarginalizationMethod(BCMarginalizationMethod type);
780 
784  std::string DumpCurrentMarginalizationMethod()
785  { return DumpMarginalizationMethod(fMarginalizationMethodCurrent); }
786 
790  std::string DumpUsedMarginalizationMethod()
791  { return DumpMarginalizationMethod(fMarginalizationMethodUsed); }
792 
797  std::string DumpOptimizationMethod(BCOptimizationMethod type);
798 
802  std::string DumpCurrentOptimizationMethod()
803  { return DumpOptimizationMethod(fOptimizationMethodCurrent); }
804 
808  std::string DumpUsedOptimizationMethod()
809  { return DumpOptimizationMethod(fOptimizationMethodUsed); }
810 
815  std::string DumpCubaIntegrationMethod(BCCubaMethod type);
816 
820  std::string DumpCubaIntegrationMethod()
821  { return DumpCubaIntegrationMethod(fCubaIntegrationMethod); }
822 
825  void SetBestFitParameters(const std::vector<double> &x)
826  { fBestFitParameters = x; }
827 
832  void SetBestFitParameters(const std::vector<double> &x, const double &new_value, double &old_value);
833 
837  unsigned GetNIntegrationVariables();
838 
842  double CalculateIntegrationVolume();
843 
846  bool CheckMarginalizationAvailability(BCMarginalizationMethod type);
847 
850  bool CheckMarginalizationIndices(TH1* hist, const std::vector<unsigned> &index);
851 
854 protected:
855 
858  int fID;
859 
862  TMinuit * fMinuit;
863 
864  double fMinuitArglist[2];
865  int fMinuitErrorFlag;
866 
869  bool fFlagIgnorePrevOptimization;
870 
873  double fSAT0;
874 
877  double fSATmin;
878 
881  TTree * fTreeSA;
882 
885  bool fFlagWriteSAToFile;
886 
887  int fSANIterations;
888  double fSATemperature;
889  double fSALogProb;
890  std::vector<double> fSAx;
891 
894  std::vector<BCH1D*> fMarginalized1D;
895 
898  std::vector<BCH2D*> fMarginalized2D;
899 
900  protected:
903  unsigned IntegrationOutputFrequency() const;
904 
910  void LogOutputAtStartOfIntegration(BCIntegrationMethod type, BCCubaMethod cubatype);
911  void LogOutputAtIntegrationStatusUpdate(BCIntegrationMethod type, double integral, double absprecision, int nIterations);
912  void LogOutputAtEndOfIntegration(double integral, double absprecision, double relprecision, int nIterations);
913 
917  void Copy(const BCIntegrate & bcintegrate);
918 
921  bool fFlagMarginalized;
922 
923 private:
924 
935  std::vector<double> FindModeMinuit(std::vector<double> &mode, std::vector<double> &errors, std::vector<double> start = std::vector<double>(0), int printlevel = 1);
936 
944  std::vector<double> FindModeMCMC(std::vector<double> &mode, std::vector<double> &errors);
945 
956  std::vector<double> FindModeSA(std::vector<double> &mode, std::vector<double> &errors, std::vector<double> start = std::vector<double>(0));
957 
961  double IntegrateCuba()
962  { return IntegrateCuba(fCubaIntegrationMethod); }
963 
968  double IntegrateCuba(BCCubaMethod cubatype);
969 
973  double IntegrateSlice();
974 
977  BCIntegrate::BCOptimizationMethod fOptimizationMethodCurrent;
978 
982  BCIntegrate::BCOptimizationMethod fOptimizationMethodUsed;
983 
986  BCIntegrate::BCIntegrationMethod fIntegrationMethodCurrent;
987 
990  BCIntegrate::BCIntegrationMethod fIntegrationMethodUsed;
991 
994  BCIntegrate::BCMarginalizationMethod fMarginalizationMethodCurrent;
995 
998  BCIntegrate::BCMarginalizationMethod fMarginalizationMethodUsed;
999 
1002  BCIntegrate::BCSASchedule fSASchedule;
1003 
1006  unsigned fNIterationsMin;
1007 
1010  unsigned fNIterationsMax;
1011 
1014  unsigned fNIterationsPrecisionCheck;
1015 
1018  unsigned fNIterationsOutput;
1019 
1022  int fNIterations;
1023 
1026  std::vector<double> fBestFitParameters;
1027 
1030  std::vector<double> fBestFitParameterErrors;
1031 
1034  double fLogMaximum;
1035 
1038  double fIntegral;
1039 
1041  double fRelativePrecision;
1042 
1044  double fAbsolutePrecision;
1045 
1048  double fError;
1049 
1051  BCCubaMethod fCubaIntegrationMethod;
1052  BCCubaOptions::Vegas fCubaVegasOptions;
1053  BCCubaOptions::Suave fCubaSuaveOptions;
1054  BCCubaOptions::Divonne fCubaDivonneOptions;
1055  BCCubaOptions::Cuhre fCubaCuhreOptions;
1056 
1057 };
1058 
1059 // ---------------------------------------------------------
1060 
1061 #endif
BCIntegrate
A class for handling numerical operations for models.
BCH2D
A class for handling 2D distributions.
Definition: BCH2D.h:37
BCParameter
A class representing a parameter of a model.
Definition: BCParameter.h:28
BCEngineMCMC
An engine class for Markov Chain Monte Carlo.
Definition: BCEngineMCMC.h:41
BCH1D
A class for handling 1D distributions.
Definition: BCH1D.h:31