00001 #ifndef __BCINTEGRATE__H 00002 #define __BCINTEGRATE__H 00003 00004 /*! 00005 * \class BCIntegrate 00006 * \brief A class for handling numerical operations for models. 00007 * \author Daniel Kollar 00008 * \author Kevin Kröninger 00009 * \version 1.0 00010 * \date 08.2008 00011 * \detail This is a base class for a model class. It contains 00012 * numerical methods to carry out the integration, marginalization, 00013 * peak finding etc. 00014 */ 00015 00016 /* 00017 * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger. 00018 * All rights reserved. 00019 * 00020 * For the licensing terms see doc/COPYING. 00021 */ 00022 00023 // --------------------------------------------------------- 00024 00025 #include <vector> 00026 00027 #include <math.h> 00028 00029 #include "BAT/BCEngineMCMC.h" 00030 #include "BAT/BCParameter.h" 00031 #include "BAT/BCDataPoint.h" 00032 00033 // ROOT classes 00034 class TH1D; 00035 class TH2D; 00036 class TRandom3; 00037 class TMinuit; 00038 class TTree; 00039 00040 00041 #define DEBUG 0 00042 00043 // --------------------------------------------------------- 00044 00045 class BCIntegrate : public BCEngineMCMC 00046 { 00047 00048 public: 00049 00050 /** \name Enumerators */ 00051 /* @{ */ 00052 00053 /** 00054 * An enumerator for the integration algorithm */ 00055 enum BCIntegrationMethod { kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba }; 00056 00057 /** 00058 * An enumerator for the marginalization algorithm */ 00059 enum BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis }; 00060 00061 /** 00062 * An enumerator for the mode finding algorithm */ 00063 enum BCOptimizationMethod { kOptSA, kOptMetropolis, kOptMinuit }; 00064 00065 /** 00066 * An enumerator for the Simulated Annealing schedule */ 00067 enum BCSASchedule { kSACauchy, kSABoltzmann, kSACustom }; 00068 00069 /* @} */ 00070 00071 /** \name Constructors and destructors */ 00072 /* @{ */ 00073 00074 /** 00075 * The default constructor */ 00076 BCIntegrate(); 00077 00078 /** 00079 * A constructor */ 00080 BCIntegrate(BCParameterSet * par); 00081 00082 /** 00083 * The default destructor */ 00084 virtual ~BCIntegrate(); 00085 00086 /* @} */ 00087 00088 /** \name Member functions (get) */ 00089 /* @{ */ 00090 00091 /** 00092 * @return The integration method */ 00093 BCIntegrate::BCIntegrationMethod GetIntegrationMethod() 00094 { return fIntegrationMethod; }; 00095 00096 /** 00097 * @return The marginalization method */ 00098 BCIntegrate::BCMarginalizationMethod GetMarginalizationMethod() 00099 { return fMarginalizationMethod; }; 00100 00101 /** 00102 * @return The current optimization method */ 00103 BCIntegrate::BCOptimizationMethod GetOptimizationMethod() 00104 { return fOptimizationMethod; }; 00105 00106 /** 00107 * @return The optimization method used to find the mode */ 00108 BCIntegrate::BCOptimizationMethod GetOptimizationMethodMode() 00109 { return fOptimizationMethodMode; }; 00110 00111 /** 00112 * @return The Simulated Annealing schedule */ 00113 BCIntegrate::BCSASchedule GetSASchedule() 00114 { return fSASchedule; }; 00115 00116 /** 00117 * Fills a vector of random numbers between 0 and 1 into a vector 00118 * @param A vector of doubles */ 00119 void GetRandomVector(std::vector <double> &x); 00120 00121 virtual void GetRandomVectorMetro(std::vector <double> &x); 00122 00123 /** 00124 * Fills a vector of (flat) random numbers in the limits of the parameters and returns 00125 * the probability at that point 00126 * @param x A vector of doubles 00127 * @return The (unnormalized) probability at the random point */ 00128 double GetRandomPoint(std::vector <double> &x); 00129 00130 /** 00131 * Fills a vector of random numbers in the limits of the parameters sampled by the sampling 00132 * function and returns the probability at that point 00133 * @param x A vector of doubles 00134 * @return The (unnormalized) probability at the random point */ 00135 double GetRandomPointImportance(std::vector <double> &x); 00136 00137 /** 00138 * Fills a vector of random numbers in the limits of the parameters sampled by the probality 00139 * function and returns the probability at that point (Metropolis) 00140 * @param x A vector of doubles */ 00141 void GetRandomPointMetro(std::vector <double> &x); 00142 00143 /** 00144 * Fills a vector of random numbers in the limits of the parameters sampled by the sampling 00145 * function and returns the probability at that point (Metropolis) 00146 * @param x A vector of doubles */ 00147 void GetRandomPointSamplingMetro(std::vector <double> &x); 00148 00149 /** 00150 * @return The number of iterations per dimension for the Monte Carlo integration */ 00151 int GetNiterationsPerDimension() 00152 { return fNiterPerDimension; }; 00153 00154 /** 00155 * @return Number of samples per 2D bin per variable in the Metropolis marginalization. */ 00156 int GetNSamplesPer2DBin() 00157 { return fNSamplesPer2DBin; }; 00158 00159 /** 00160 * @return The number of variables to integrate over */ 00161 int GetNvar() 00162 { return fNvar; }; 00163 00164 /** 00165 * @return The number of maximum iterations for Monte Carlo integration */ 00166 int GetNIterationsMax() 00167 { return fNIterationsMax; }; 00168 00169 /** 00170 * @return The number of iterations for the most recent Monte Carlo integration */ 00171 int GetNIterations() 00172 { return fNIterations; }; 00173 00174 /** 00175 * @return The relative precision for numerical integration */ 00176 double GetRelativePrecision() 00177 { return fRelativePrecision; }; 00178 00179 /* 00180 * @return The uncertainty in the most recent Monte Carlo integration */ 00181 double GetError() 00182 { return fError; }; 00183 00184 /* 00185 * @return number of bins per dimension for the marginalized distributions */ 00186 int GetNbins() 00187 { return fNbins; }; 00188 00189 /* 00190 * @return Minuit used for mode finding */ 00191 TMinuit * GetMinuit(); 00192 00193 int GetMinuitErrorFlag() 00194 { return fMinuitErrorFlag; }; 00195 00196 /* 00197 * @return The ROOT tree containing the Markov chain */ 00198 TTree * GetMarkovChainTree() 00199 { return fMarkovChainTree; }; 00200 00201 /* 00202 * Returns the actual point in the markov chain */ 00203 std::vector<double> * GetMarkovChainPoint() 00204 { return &fXmetro1; }; 00205 00206 /** 00207 * Returns the iteration of the MCMC */ 00208 int * GetMCMCIteration() 00209 { return &fMCMCIteration; }; 00210 00211 /** 00212 * Returns the value of the loglikelihood at the point fXmetro1 */ 00213 double * GetMarkovChainValue() 00214 { return &fMarkovChainValue; }; 00215 00216 /** 00217 * Returns the Simulated Annealing starting temperature. */ 00218 double GetSAT0() 00219 { return fSAT0; } 00220 00221 /** 00222 * Returns the Simulated Annealing threshhold temperature. */ 00223 double GetSATmin() 00224 { return fSATmin; } 00225 00226 double * GetSAP2Temperature() 00227 { return &fSATemperature; }; 00228 00229 double * GetSAP2LogProb() 00230 { return &fSALogProb; }; 00231 00232 std::vector<double> * GetSAP2x() 00233 { return &fSAx; }; 00234 00235 int * GetSAP2NIterations() 00236 { return &fSANIterations; }; 00237 00238 /* @} */ 00239 00240 /** \name Member functions (set) */ 00241 /* @{ */ 00242 00243 void SetMinuitArlist(double * arglist) 00244 { fMinuitArglist[0] = arglist[0]; 00245 fMinuitArglist[1] = arglist[1]; }; 00246 00247 void SetFlagIgnorePrevOptimization(bool flag) 00248 { fFlagIgnorePrevOptimization = flag; }; 00249 00250 /** 00251 * @param par The parameter set which gets translated into array 00252 * needed for the Monte Carlo integration */ 00253 void SetParameters(BCParameterSet * par); 00254 00255 /** 00256 * @param varlist A list of parameters */ 00257 void SetVarList(int * varlist); 00258 00259 /** 00260 * @param index The index of the variable to be set */ 00261 void SetVar(int index){fVarlist[index]=1;}; 00262 00263 /** 00264 * @param method The integration method */ 00265 void SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method); 00266 00267 /** 00268 * @param method The marginalization method */ 00269 void SetMarginalizationMethod(BCIntegrate::BCMarginalizationMethod method) 00270 { fMarginalizationMethod = method; }; 00271 00272 /** 00273 * @param method The mode finding method */ 00274 void SetOptimizationMethod(BCIntegrate::BCOptimizationMethod method) 00275 { fOptimizationMethod = method; }; 00276 00277 /** 00278 * @param method The Simulated Annealing schedule */ 00279 void SetSASchedule(BCIntegrate::BCSASchedule schedule) 00280 { fSASchedule = schedule; }; 00281 00282 /** 00283 * @param niterations Number of iterations per dimension for Monte Carlo integration. */ 00284 void SetNiterationsPerDimension(int niterations) 00285 { fNiterPerDimension = niterations; }; 00286 00287 /** 00288 * @param n Number of samples per 2D bin per variable in the Metropolis marginalization. 00289 * Default is 100. */ 00290 void SetNSamplesPer2DBin(int n) 00291 { fNSamplesPer2DBin = n; }; 00292 00293 /** 00294 * @param niterations The maximum number of iterations for Monte Carlo integration */ 00295 void SetNIterationsMax(int niterations) 00296 { fNIterationsMax = niterations; }; 00297 00298 /** 00299 * @param relprecision The relative precision envisioned for Monte 00300 * Carlo integration */ 00301 void SetRelativePrecision(double relprecision) 00302 { fRelativePrecision = relprecision; }; 00303 00304 /** 00305 * @param n Number of bins per dimension for the marginalized 00306 * distributions. Default is 100. Minimum number allowed is 2. */ 00307 00308 /** 00309 * Set the number of bins for the marginalized distribution of a 00310 * parameter. 00311 * @param nbins Number of bins (default = 100) 00312 * @param index Index of the parameter. */ 00313 void SetNbins(int nbins, int index = -1); 00314 00315 // void SetNbins(int n); 00316 // void SetNbinsX(int n); 00317 // void SetNbinsY(int n); 00318 00319 /* 00320 * Turn on or off the filling of the error band during the MCMC run. 00321 * @param flag set to true for turning on the filling */ 00322 void SetFillErrorBand(bool flag = true) 00323 { fFillErrorBand=flag; }; 00324 00325 /* 00326 * Turn off filling of the error band during the MCMC run. 00327 * This method is equivalent to SetFillErrorBand(false) */ 00328 void UnsetFillErrorBand() 00329 { fFillErrorBand=false; }; 00330 00331 /** 00332 * Sets index of the x values in function fits. 00333 * @param index Index of the x values */ 00334 void SetFitFunctionIndexX(int index) 00335 { fFitFunctionIndexX = index; }; 00336 00337 /** 00338 * Sets index of the y values in function fits. @param index Index 00339 * of the y values */ 00340 void SetFitFunctionIndexY(int index) 00341 { fFitFunctionIndexY = index; }; 00342 00343 void SetFitFunctionIndices(int indexx, int indexy) 00344 { this -> SetFitFunctionIndexX(indexx); 00345 this -> SetFitFunctionIndexY(indexy); }; 00346 00347 /** 00348 * Sets the data point containing the lower boundaries of possible 00349 * data values */ 00350 void SetDataPointLowerBoundaries(BCDataPoint * datasetlowerboundaries) 00351 { fDataPointLowerBoundaries = datasetlowerboundaries; }; 00352 00353 /** 00354 * Sets the data point containing the upper boundaries of possible 00355 * data values */ 00356 void SetDataPointUpperBoundaries(BCDataPoint * datasetupperboundaries) 00357 { fDataPointUpperBoundaries = datasetupperboundaries; }; 00358 00359 /** 00360 * Sets the lower boundary of possible data values for a particular 00361 * variable */ 00362 void SetDataPointLowerBoundary(int index, double lowerboundary) 00363 { fDataPointLowerBoundaries -> SetValue(index, lowerboundary); }; 00364 00365 /** 00366 * Sets the upper boundary of possible data values for a particular 00367 * variable */ 00368 void SetDataPointUpperBoundary(int index, double upperboundary) 00369 { fDataPointUpperBoundaries -> SetValue(index, upperboundary); }; 00370 00371 /** 00372 * Flag for writing Markov chain to ROOT file (true) or not (false) 00373 */ 00374 void WriteMarkovChain(bool flag) 00375 { fFlagWriteMarkovChain = flag; fMCMCFlagWriteChainToFile = flag;}; 00376 00377 /** 00378 * Sets the ROOT tree containing the Markov chain */ 00379 void SetMarkovChainTree(TTree * tree) 00380 { fMarkovChainTree = tree; }; 00381 00382 /** 00383 * Sets the initial position for the Markov chain */ 00384 void SetMarkovChainInitialPosition(std::vector<double> position); 00385 00386 /** 00387 * Sets the step size for Markov chains */ 00388 void SetMarkovChainStepSize(double stepsize) 00389 { fMarkovChainStepSize = stepsize; }; 00390 00391 /** 00392 * Sets the number of iterations in the markov chain */ 00393 void SetMarkovChainNIterations(int niterations) 00394 { fMarkovChainNIterations = niterations; 00395 fMarkovChainAutoN = false; } 00396 00397 /** 00398 * Sets a flag for automatically calculating the number of iterations */ 00399 void SetMarkovChainAutoN(bool flag) 00400 { fMarkovChainAutoN = flag; }; 00401 00402 /** 00403 * Sets mode */ 00404 void SetMode(std::vector <double> mode); 00405 00406 /** 00407 * Sets errorband histogram */ 00408 void SetErrorBandHisto(TH2D * h) 00409 { fErrorBandXY = h; }; 00410 00411 /** 00412 * @param T0 new value for Simulated Annealing starting temperature. */ 00413 void SetSAT0(double T0) 00414 { fSAT0 = T0; } 00415 00416 /** 00417 * @param Tmin new value for Simulated Annealing threshold temperature. */ 00418 void SetSATmin(double Tmin) 00419 { fSATmin = Tmin; } 00420 00421 void SetFlagWriteSAToFile(bool flag) 00422 { fFlagWriteSAToFile = flag; }; 00423 00424 /* 00425 * Sets the tree containing the Markov chains. */ 00426 void SetSATree(TTree* tree) 00427 { fTreeSA = tree; }; 00428 00429 /* @} */ 00430 00431 /** \name Member functions (miscellaneous methods) */ 00432 /* @{ */ 00433 00434 /** 00435 * Frees the memory for integration variables */ 00436 void DeleteVarList(); 00437 00438 /** 00439 * Sets all values of the variable list to a particular value 00440 * @v The value */ 00441 void ResetVarlist(int v); 00442 00443 /** 00444 * Set value of a particular integration variable to 0. 00445 * @param index The index of the variable */ 00446 void UnsetVar(int index) 00447 { fVarlist[index] = 0; }; 00448 00449 /** 00450 * Evaluate the unnormalized probability at a point in parameter space. 00451 * Method needs to be overloaded by the user. 00452 * @param x The point in parameter space 00453 * @return The unnormalized probability */ 00454 virtual double Eval(std::vector <double> x); 00455 00456 /** 00457 * Evaluate the natural logarithm of the Eval function. For better numerical 00458 * stability, this method should also be overloaded by the user. 00459 * @param x The point in parameter space 00460 * @return log(Eval(x)) */ 00461 virtual double LogEval(std::vector <double> x); 00462 00463 /** 00464 * Evaluate the sampling function at a point in parameter space. 00465 * Method needs to be overloaded by the user. 00466 * @param x The point in parameter space 00467 * @return The value of the sampling function */ 00468 virtual double EvalSampling(std::vector <double> x); 00469 00470 /** 00471 * Evaluate the natural logarithm of the EvalSampling function. 00472 * Method needs to be overloaded by the user. 00473 * @param x The point in parameter space 00474 * @return log(EvalSampling(x)) */ 00475 double LogEvalSampling(std::vector <double> x); 00476 00477 /** 00478 * Evaluate the un-normalized probability at a point in parameter space 00479 * and prints the result to the log. 00480 * @param x The point in parameter space 00481 * @return The un-normalized probability 00482 * @see Eval(std::vector <double> x) */ 00483 double EvalPrint(std::vector <double> x); 00484 00485 /** 00486 * Defines a fit function. 00487 * @param parameters A set of parameter values 00488 * @param x A vector of x-values 00489 * @return The value of the fit function at the x-values given a set of parameters */ 00490 virtual double FitFunction(std::vector <double> x, std::vector <double> parameters) 00491 { return 0.; }; 00492 00493 /** 00494 * Does the integration over the un-normalized probability. 00495 * @return The normalization value */ 00496 double Integrate(); 00497 00498 /** 00499 * Perfoms the Monte Carlo integration. For details see documentation. 00500 * @param x An initial point in parameter space 00501 * @param varlist A list of variables 00502 * @return The integral */ 00503 double IntegralMC(std::vector <double> x, int * varlist); 00504 00505 double IntegralMC(std::vector <double> x); 00506 00507 /** 00508 * Perfoms the Metropolis Monte Carlo integration. For details see documentation. 00509 * @param x An initial point in parameter space 00510 * @return The integral */ 00511 double IntegralMetro(std::vector <double> x); 00512 00513 /** 00514 * Perfoms the importance sampling Monte Carlo integration. For details see documentation. 00515 * @param x An initial point in parameter space 00516 * @return The integral */ 00517 double IntegralImportance(std::vector <double> x); 00518 00519 /** 00520 * Calculate integral using the Cuba library. For details see documentation. 00521 * @param method A short cut for the method 00522 * @param parameters_double A vector of parameters (double) 00523 * @param parameters_int A vector of parameters (int) 00524 * @return The integral */ 00525 double CubaIntegrate(int method, std::vector<double> parameters_double, std::vector<double> parameters_int); 00526 00527 double CubaIntegrate(); 00528 00529 /** 00530 * Integrand for the Cuba library. 00531 * @param ndim The number of dimensions to integrate over 00532 * @param xx The point in parameter space to integrate over (scaled to 0 - 1 per dimension) 00533 * @param ncomp The number of components of the integrand (usually 1) 00534 * @param ff The function value 00535 * @return The integral */ 00536 static void CubaIntegrand(const int * ndim, const double xx[], const int * ncomp, double ff[]); 00537 00538 /** 00539 * Performs the marginalization with respect to one parameter. 00540 * @param parameter The parameter w.r.t. which the marginalization is performed 00541 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00542 TH1D * Marginalize(BCParameter * parameter); 00543 00544 /** 00545 * Performs the marginalization with respect to two parameters. 00546 * @param parameter1 The first parameter w.r.t. which the marginalization is performed 00547 * @param parameter2 The second parameter w.r.t. which the marginalization is performed 00548 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00549 TH2D * Marginalize(BCParameter * parameter1, BCParameter * parameter2); 00550 00551 /** 00552 * Performs the marginalization with respect to one parameter using 00553 * the simple Monte Carlo technique. 00554 * @param parameter The parameter w.r.t. which the marginalization is performed 00555 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00556 TH1D * MarginalizeByIntegrate(BCParameter * parameter); 00557 00558 /** 00559 * Performs the marginalization with respect to two parameters using 00560 * the simple Monte Carlo technique. 00561 * @param parameter1 The first parameter w.r.t. which the marginalization is performed 00562 * @param parameter2 The second parameter w.r.t. which the marginalization is performed 00563 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00564 TH2D * MarginalizeByIntegrate(BCParameter * parameter1, BCParameter * parameter2); 00565 00566 /** 00567 * Performs the marginalization with respect to one parameter using 00568 * the Metropolis algorithm. 00569 * @param parameter The parameter w.r.t. which the marginalization is performed 00570 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00571 TH1D * MarginalizeByMetro(BCParameter * parameter); 00572 00573 /** 00574 * Performs the marginalization with respect to two parameters using the Metropolis algorithm. 00575 * @param parameter1 The first parameter w.r.t. which the marginalization is performed 00576 * @param parameter2 The second parameter w.r.t. which the marginalization is performed 00577 * @return A histogram which contains the marginalized probability distribution (normalized to 1) */ 00578 TH2D * MarginalizeByMetro(BCParameter * parameter1, BCParameter * parameter2); 00579 00580 /** 00581 * Performs the marginalization with respect to every single parameter as well as with respect 00582 * all combinations to two parameters using the Metropolis algorithm. 00583 * @param name Basename for the histograms (e.g. model name) 00584 * @return Total number of marginalized distributions */ 00585 int MarginalizeAllByMetro(const char * name); 00586 00587 /** 00588 * @param parIndex1 Index of parameter 00589 * @return Pointer to 1D histogram (TH1D) of marginalized distribution wrt. parameter with given index. 00590 */ 00591 TH1D * GetH1D(int parIndex); 00592 00593 /** 00594 * @param parIndex1 Index of first parameter 00595 * @param parIndex2 Index of second parameter, with parIndex2>parIndex1 00596 * @return Index of the distribution in the vector of 2D distributions, which corresponds 00597 * to the combination of parameters with given indeces */ 00598 int GetH2DIndex(int parIndex1, int parIndex2); 00599 00600 /** 00601 * @param parIndex1 Index of first parameter 00602 * @param parIndex2 Index of second parameter, with parIndex2>parIndex1 00603 * @return Pointer to 2D histogram (TH2D) of marginalized distribution wrt. parameters with given indeces. 00604 */ 00605 TH2D * GetH2D(int parIndex1, int parIndex2); 00606 00607 /** 00608 * Initializes the Metropolis algorithm (for details see manual) */ 00609 void InitMetro(); 00610 00611 void SAInitialize(); 00612 00613 /** 00614 * Does the mode finding */ 00615 // void FindMode(); 00616 00617 /** 00618 * Does the mode finding using Minuit. If starting point is not specified, 00619 * finding will start from the center of the parameter space. 00620 * @param start point in parameter space from which the mode finding is started. 00621 * @param printlevel The print level. */ 00622 virtual void FindModeMinuit(std::vector<double> start = std::vector<double>(0), int printlevel = 1); 00623 00624 /** 00625 * Does the mode finding using Markov Chain Monte Carlo */ 00626 void FindModeMCMC(); 00627 // void FindModeMCMC(int flag_run = 0); 00628 00629 /** 00630 * Does the mode finding using Simulated Annealing. If starting point 00631 * is not specified, finding will start from the center of the 00632 * parameter space. 00633 * @param start point in parameter space from thich the mode finding is started. */ 00634 void FindModeSA(std::vector<double> start = std::vector<double>(0)); 00635 00636 /** 00637 * Temperature annealing schedule for use with Simulated Annealing. 00638 * Delegates to the appropriate method according to 00639 * fSASchedule. 00640 * @param t iterator for lowering the temperature over time. */ 00641 double SATemperature(double t); 00642 00643 /** 00644 * Temperature annealing schedule for use with Simulated Annealing. 00645 * This method is used for Boltzmann annealing schedule. 00646 * @param t iterator for lowering the temperature over time. */ 00647 double SATemperatureBoltzmann(double t); 00648 00649 /** 00650 * Temperature annealing schedule for use with Simulated Annealing. 00651 * This method is used for Cauchy annealing schedule. 00652 * @param t iterator for lowering the temperature over time. */ 00653 double SATemperatureCauchy(double t); 00654 00655 /** 00656 * Temperature annealing schedule for use with Simulated Annealing. 00657 * This is a virtual method to be overridden by a user-defined 00658 * custom temperature schedule. 00659 * @param t iterator for lowering the temperature over time. */ 00660 virtual double SATemperatureCustom(double t); 00661 00662 /** 00663 * Generates a new state in a neighbourhood around x that is to be 00664 * accepted or rejected by the Simulated Annealing algorithm. 00665 * Delegates the generation to the appropriate method according 00666 * to fSASchedule. 00667 * @param x last state. 00668 * @param t time iterator to determine current temperature. */ 00669 std::vector<double> GetProposalPointSA(std::vector<double> x, int t); 00670 00671 /** 00672 * Generates a new state in a neighbourhood around x that is to be 00673 * accepted or rejected by the Simulated Annealing algorithm. 00674 * This method is used for Boltzmann annealing schedule. 00675 * @param x last state. 00676 * @param t time iterator to determine current temperature. */ 00677 std::vector<double> GetProposalPointSABoltzmann(std::vector<double> x, int t); 00678 00679 /** 00680 * Generates a new state in a neighbourhood around x that is to be 00681 * accepted or rejected by the Simulated Annealing algorithm. 00682 * This method is used for Cauchy annealing schedule. 00683 * @param x last state. 00684 * @param t time iterator to determine current temperature. */ 00685 std::vector<double> GetProposalPointSACauchy(std::vector<double> x, int t); 00686 00687 /** 00688 * Generates a new state in a neighbourhood around x that is to be 00689 * accepted or rejected by the Simulated Annealing algorithm. 00690 * This is a virtual method to be overridden by a user-defined 00691 * custom proposal function. 00692 * @param x last state. 00693 * @param t time iterator to determine current temperature. */ 00694 virtual std::vector<double> GetProposalPointSACustom(std::vector<double> x, int t); 00695 00696 /** 00697 * Generates a uniform distributed random point on the surface of 00698 * a fNvar-dimensional Hypersphere. 00699 * Used as a helper to generate proposal points for Cauchy annealing. */ 00700 std::vector<double> SAHelperGetRandomPointOnHypersphere(); 00701 00702 /** 00703 * Generates the radial part of a n-dimensional Cauchy distribution. 00704 * Helper function for Cauchy annealing. */ 00705 double SAHelperGetRadialCauchy(); 00706 00707 /** 00708 * Returns the Integral of sin^dim from 0 to theta. 00709 * Helper function needed for generating Cauchy distributions. */ 00710 double SAHelperSinusToNIntegral(int dim, double theta); 00711 00712 00713 static void FCNLikelihood(int &npar, double * grad, double &fval, double * par, int flag); 00714 00715 /* 00716 * Method executed for every iteration of the MCMC. User's code should be 00717 * provided via overloading in the derived class*/ 00718 virtual void MCMCUserIterationInterface() 00719 {}; 00720 00721 /* @} */ 00722 00723 private: 00724 00725 /** 00726 * Set of parameters for the integration. */ 00727 BCParameterSet * fx; 00728 00729 /** 00730 * Array containing the lower boundaries of the variables to integrate over. */ 00731 double * fMin; 00732 00733 /** 00734 * Array containing the upper boundaries of the variables to integrate over. */ 00735 double * fMax; 00736 00737 /** 00738 * List of variables containing a flag whether to integrate over them or not. */ 00739 int * fVarlist; 00740 00741 /** 00742 * Number of iteration per dimension for Monte Carlo integration. */ 00743 int fNiterPerDimension; 00744 00745 /** 00746 * Current integration method */ 00747 BCIntegrate::BCIntegrationMethod fIntegrationMethod; 00748 00749 /** 00750 * Current marginalization method */ 00751 BCIntegrate::BCMarginalizationMethod fMarginalizationMethod; 00752 00753 /** 00754 * Current mode finding method */ 00755 BCIntegrate::BCOptimizationMethod fOptimizationMethod; 00756 00757 /** 00758 * Method with which the global mode was found (can differ from 00759 * fOptimization method in case more than one algorithm is used). */ 00760 BCIntegrate::BCOptimizationMethod fOptimizationMethodMode; 00761 00762 /** 00763 * Current Simulated Annealing schedule */ 00764 BCIntegrate::BCSASchedule fSASchedule; 00765 00766 /** 00767 * Maximum number of iterations */ 00768 int fNIterationsMax; 00769 00770 /** 00771 * Number of iterations in the most recent Monte Carlo integation */ 00772 int fNIterations; 00773 00774 /** 00775 * Relative precision aimed at in the Monte Carlo integation */ 00776 double fRelativePrecision; 00777 00778 /** 00779 * The uncertainty in the most recent Monte Carlo integration */ 00780 double fError; 00781 00782 /** 00783 * The number of iterations in the Metropolis integration */ 00784 int fNmetro; 00785 int fNacceptedMCMC; 00786 00787 /** 00788 * A vector of points in parameter space used for the Metropolis algorithm */ 00789 std::vector <double> fXmetro0; 00790 00791 /** 00792 * A vector of points in parameter space used for the Metropolis algorithm */ 00793 std::vector <double> fXmetro1; 00794 00795 /** 00796 * A double containing the log likelihood value at the point fXmetro1 */ 00797 double fMarkovChainValue; 00798 00799 /* 00800 * Method executed for every iteration of the MCMC, overloaded from BCEngineMCMC. */ 00801 void MCMCIterationInterface(); 00802 00803 /* 00804 * Fill error band histogram for curreent iteration. This method is called from MCMCIterationInterface() */ 00805 void MCMCFillErrorBand(); 00806 00807 protected: 00808 00809 /** 00810 * Number of variables to integrate over. */ 00811 int fNvar; 00812 00813 /** 00814 * Number of bins per dimension for the marginalized distributions */ 00815 int fNbins; 00816 00817 /** 00818 * Number of samples per 2D bin per variable in the Metropolis 00819 * marginalization. */ 00820 int fNSamplesPer2DBin; 00821 00822 /** 00823 * Step size in the Markov chain relative to min and max */ 00824 double fMarkovChainStepSize; 00825 00826 int fMarkovChainNIterations; 00827 00828 bool fMarkovChainAutoN; 00829 00830 /** 00831 * data point containing the lower boundaries of possible data values */ 00832 BCDataPoint * fDataPointLowerBoundaries; 00833 00834 /** 00835 * data point containing the upper boundaries of possible data values */ 00836 BCDataPoint * fDataPointUpperBoundaries; 00837 00838 std::vector <bool> fDataFixedValues; 00839 00840 /** 00841 * A ROOT random number generator */ 00842 TRandom3 * fRandom; 00843 00844 /** 00845 * A vector of best fit parameters estimated from the global 00846 * probability and the estimate on their uncertainties */ 00847 std::vector <double> fBestFitParameters; 00848 std::vector <double> fBestFitParameterErrors; 00849 00850 /** 00851 * A vector of best fit parameters estimated from the marginalized probability */ 00852 std::vector <double> fBestFitParametersMarginalized; 00853 00854 /** 00855 * Vector of TH1D histograms for marginalized probability distributions */ 00856 std::vector <TH1D *> fHProb1D; 00857 00858 /** 00859 * Vector of TH2D histograms for marginalized probability distributions */ 00860 std::vector <TH2D *> fHProb2D; 00861 00862 /** 00863 * Flag whether or not to fill the error band */ 00864 bool fFillErrorBand; 00865 00866 /** 00867 * The indeces for function fits */ 00868 int fFitFunctionIndexX; 00869 int fFitFunctionIndexY; 00870 00871 /** 00872 * A flag for single point evaluation of the error "band" */ 00873 bool fErrorBandContinuous; 00874 std::vector <double> fErrorBandX; 00875 00876 /** 00877 * The error band histogram */ 00878 TH2D * fErrorBandXY; 00879 00880 /** 00881 * Number of X bins of the error band histogram */ 00882 int fErrorBandNbinsX; 00883 00884 /** 00885 * Nnumber of Y bins of the error band histogram */ 00886 int fErrorBandNbinsY; 00887 00888 /** 00889 * Minuit */ 00890 TMinuit * fMinuit; 00891 00892 double fMinuitArglist[2]; 00893 int fMinuitErrorFlag; 00894 00895 /** 00896 * Flag for ignoring older results of minimization */ 00897 double fFlagIgnorePrevOptimization; 00898 00899 /** 00900 * Flag for writing Markov chain to file */ 00901 bool fFlagWriteMarkovChain; 00902 00903 /** 00904 * ROOT tree containing the Markov chain */ 00905 TTree * fMarkovChainTree; 00906 00907 /** 00908 * Iteration of the MCMC */ 00909 int fMCMCIteration; 00910 00911 /** 00912 * Starting temperature for Simulated Annealing */ 00913 double fSAT0; 00914 00915 /** 00916 * Minimal/Threshold temperature for Simulated Annealing */ 00917 double fSATmin; 00918 00919 /** 00920 * Tree for the Simulated Annealing */ 00921 TTree * fTreeSA; 00922 00923 /** 00924 * Flag deciding whether to write SA to file or not. */ 00925 bool fFlagWriteSAToFile; 00926 00927 int fSANIterations; 00928 double fSATemperature; 00929 double fSALogProb; 00930 std::vector<double> fSAx; 00931 00932 }; 00933 00934 // --------------------------------------------------------- 00935 00936 #endif