Reference Manual

Go to the documentation of this file.
00001 #include "../interface/CascadeMinimizer.h"
00002 #include "../interface/ProfiledLikelihoodRatioTestStatExt.h"
00003 #include "../interface/ProfileLikelihood.h"
00004 #include "../interface/CloseCoutSentry.h"
00005 #include "../interface/utils.h"
00006 
00007 #include <Math/MinimizerOptions.h>
00008 #include <RooCategory.h>
00009 #include <RooNumIntConfig.h>
00010 
00011 
00012 boost::program_options::options_description CascadeMinimizer::options_("Cascade Minimizer options");
00013 std::vector<CascadeMinimizer::Algo> CascadeMinimizer::fallbacks_;
00014 bool CascadeMinimizer::preScan_;
00015 bool CascadeMinimizer::poiOnlyFit_;
00016 bool CascadeMinimizer::singleNuisFit_;
00017 bool CascadeMinimizer::setZeroPoint_ = true;
00018 bool CascadeMinimizer::oldFallback_ = true;
00019 
00020 CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi, int initialStrategy) :
00021     nll_(nll),
00022     minimizer_(nll_),
00023     mode_(mode),
00024     strategy_(initialStrategy),
00025     poi_(poi)
00026 {
00027 }
00028 
00029 bool CascadeMinimizer::improve(int verbose, bool cascade) 
00030 {
00031     minimizer_.setPrintLevel(verbose-2);  
00032     minimizer_.setStrategy(strategy_);
00033     bool outcome = improveOnce(verbose-2);
00034     if (cascade && !outcome && !fallbacks_.empty()) {
00035         std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType());
00036         std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo());
00037         float       nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance());
00038         int         nominalStrat(strategy_);
00039         if (verbose > 0) std::cerr << "Failed minimization with " << nominalType << "," << nominalAlgo << " and tolerance " << nominalTol << std::endl;
00040         for (std::vector<Algo>::const_iterator it = fallbacks_.begin(), ed = fallbacks_.end(); it != ed; ++it) {
00041             ProfileLikelihood::MinimizerSentry minimizerConfig(it->algo, it->tolerance != Algo::default_tolerance() ? it->tolerance : nominalTol);
00042             int myStrategy = it->strategy; if (myStrategy == Algo::default_strategy()) myStrategy = nominalStrat;
00043             if (nominalType != ROOT::Math::MinimizerOptions::DefaultMinimizerType() ||
00044                 nominalAlgo != ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo() ||
00045                 nominalTol  != ROOT::Math::MinimizerOptions::DefaultTolerance()     ||
00046                 myStrategy  != nominalStrat) {
00047                 if (verbose > 0) std::cerr << "Will fallback to minimization using " << it->algo << ", strategy " << myStrategy << " and tolerance " << it->tolerance << std::endl;
00048                 minimizer_.setStrategy(myStrategy);
00049                 outcome = improveOnce(verbose-2);
00050                 if (outcome) break;
00051             }
00052         }
00053     }
00054     if (setZeroPoint_) {
00055         cacheutils::CachingSimNLL *simnll = dynamic_cast<cacheutils::CachingSimNLL *>(&nll_);
00056         if (simnll) simnll->clearZeroPoint();
00057     }
00058     return outcome;
00059 }
00060 
00061 bool CascadeMinimizer::improveOnce(int verbose) 
00062 {
00063     std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType());
00064     std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo());
00065     if (setZeroPoint_) {
00066         cacheutils::CachingSimNLL *simnll = dynamic_cast<cacheutils::CachingSimNLL *>(&nll_);
00067         if (simnll) { 
00068             simnll->setZeroPoint();
00069         }
00070     }
00071     bool outcome = false;
00072     if (oldFallback_){
00073         outcome = nllutils::robustMinimize(nll_, minimizer_, verbose);
00074     } else {
00075         int status = minimizer_.minimize(myType.c_str(), myAlgo.c_str());
00076         outcome = (status == 0);
00077     }
00078     if (setZeroPoint_) {
00079         cacheutils::CachingSimNLL *simnll = dynamic_cast<cacheutils::CachingSimNLL *>(&nll_);
00080         if (simnll) simnll->clearZeroPoint();
00081     }
00082     return outcome;
00083 }
00084 
00085 bool CascadeMinimizer::minimize(int verbose, bool cascade) 
00086 {
00087     minimizer_.setPrintLevel(verbose-2);  
00088     minimizer_.setStrategy(strategy_);
00089     if (preScan_) minimizer_.minimize("Minuit2","Scan");
00090      // FIXME can be made smarter than this
00091     if (mode_ == Unconstrained && poiOnlyFit_) {
00092         trivialMinimize(nll_, *poi_, 200);
00093     }
00094     return improve(verbose, cascade);
00095 }
00096 
00097 
00098 void CascadeMinimizer::initOptions() 
00099 {
00100     options_.add_options()
00101         ("cminPoiOnlyFit",  "Do first a fit floating only the parameter of interest")
00102         ("cminPreScan",  "Do a scan before first minimization")
00103         ("cminSingleNuisFit", "Do first a minimization of each nuisance parameter individually")
00104         ("cminFallbackAlgo", boost::program_options::value<std::vector<std::string> >(), "Fallback algorithms if the default minimizer fails (can use multiple ones). Syntax is algo[,subalgo][,strategy][:tolerance]")
00105         ("cminSetZeroPoint", boost::program_options::value<bool>(&setZeroPoint_)->default_value(setZeroPoint_), "Change the reference point of the NLL to be zero during minimization")
00106         ("cminOldRobustMinimize", boost::program_options::value<bool>(&oldFallback_)->default_value(oldFallback_), "Use the old 'robustMinimize' logic in addition to the cascade")
00107         //("cminDefaultIntegratorEpsAbs", boost::program_options::value<double>(), "RooAbsReal::defaultIntegratorConfig()->setEpsAbs(x)")
00108         //("cminDefaultIntegratorEpsRel", boost::program_options::value<double>(), "RooAbsReal::defaultIntegratorConfig()->setEpsRel(x)")
00109         //("cminDefaultIntegrator1D", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->method1D().setLabel(x)")
00110         //("cminDefaultIntegrator1DOpen", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->method1DOpen().setLabel(x)")
00111         //("cminDefaultIntegrator2D", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->method2D().setLabel(x)")
00112         //("cminDefaultIntegrator2DOpen", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->method2DOpen().setLabel(x)")
00113         //("cminDefaultIntegratorND", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->methodND().setLabel(x)")
00114         //("cminDefaultIntegratorNDOpen", boost::program_options::value<std::string>(), "RooAbsReal::defaultIntegratorConfig()->methodNDOpen().setLabel(x)")
00115         ;
00116 }
00117 
00118 void CascadeMinimizer::applyOptions(const boost::program_options::variables_map &vm) 
00119 {
00120     using namespace std;
00121 
00122     preScan_ = vm.count("cminPreScan");
00123     poiOnlyFit_ = vm.count("cminPoiOnlyFit");
00124     singleNuisFit_ = vm.count("cminSingleNuisFit");
00125     setZeroPoint_  = vm.count("cminSetZeroPoint");
00126     if (vm.count("cminFallbackAlgo")) {
00127         vector<string> falls(vm["cminFallbackAlgo"].as<vector<string> >());
00128         for (vector<string>::const_iterator it = falls.begin(), ed = falls.end(); it != ed; ++it) {
00129             std::string algo = *it;
00130             float tolerance = Algo::default_tolerance(); 
00131             int   strategy = Algo::default_strategy(); 
00132             string::size_type idx = std::min(algo.find(";"), algo.find(":"));
00133             if (idx != string::npos && idx < algo.length()) {
00134                  tolerance = atof(algo.substr(idx+1).c_str());
00135                  algo      = algo.substr(0,idx); // DON'T SWAP THESE TWO LINES
00136             }
00137             idx = algo.find(",");
00138             if (idx != string::npos && idx < algo.length()) {
00139                 // if after the comma there's a number, then it's a strategy
00140                 if ( '0' <= algo[idx+1] && algo[idx+1] <= '9' ) {
00141                     strategy = atoi(algo.substr(idx+1).c_str());
00142                     algo     = algo.substr(0,idx); // DON'T SWAP THESE TWO LINES
00143                 } else {
00144                 // otherwise, it could be Name,subname,strategy
00145                     idx = algo.find(",",idx+1);
00146                     if (idx != string::npos && idx < algo.length()) {
00147                         strategy = atoi(algo.substr(idx+1).c_str());
00148                         algo     = algo.substr(0,idx); // DON'T SWAP THESE TWO LINES
00149                     }
00150                 }
00151             }
00152             fallbacks_.push_back(Algo(algo, tolerance, strategy));
00153             std::cout << "Configured fallback algorithm " << fallbacks_.back().algo << 
00154                             ", strategy " << fallbacks_.back().strategy   << 
00155                             ", tolerance " << fallbacks_.back().tolerance << std::endl;
00156         }
00157     }
00158     //if (vm.count("cminDefaultIntegratorEpsAbs")) RooAbsReal::defaultIntegratorConfig()->setEpsAbs(vm["cminDefaultIntegratorEpsAbs"].as<double>());
00159     //if (vm.count("cminDefaultIntegratorEpsRel")) RooAbsReal::defaultIntegratorConfig()->setEpsRel(vm["cminDefaultIntegratorEpsRel"].as<double>());
00160     //if (vm.count("cminDefaultIntegrator1D")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->method1D(), vm["cminDefaultIntegrator1D"].as<std::string>());
00161     //if (vm.count("cminDefaultIntegrator1DOpen")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->method1DOpen(), vm["cminDefaultIntegrator1DOpen"].as<std::string>());
00162     //if (vm.count("cminDefaultIntegrator2D")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->method2D(), vm["cminDefaultIntegrator2D"].as<std::string>());
00163     //if (vm.count("cminDefaultIntegrator2DOpen")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->method2DOpen(), vm["cminDefaultIntegrator2DOpen"].as<std::string>());
00164     //if (vm.count("cminDefaultIntegratorND")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->methodND(), vm["cminDefaultIntegratorND"].as<std::string>());
00165     //if (vm.count("cminDefaultIntegratorNDOpen")) setDefaultIntegrator(RooAbsReal::defaultIntegratorConfig()->methodNDOpen(), vm["cminDefaultIntegratorNDOpen"].as<std::string>());
00166 }
00167 
00168 //void CascadeMinimizer::setDefaultIntegrator(RooCategory &cat, const std::string & val) {
00169 //    if (val == "list") {
00170 //        std::cout << "States for " << cat.GetName() << std::endl;
00171 //        int i0 = cat.getBin();
00172 //        for (int i = 0, n = cat.numBins((const char *)0); i < n; ++i) {
00173 //            cat.setBin(i); std::cout << " - " << cat.getLabel() <<  ( i == i0 ? " (current default)" : "") << std::endl;
00174 //        }
00175 //        std::cout << std::endl;
00176 //        cat.setBin(i0);
00177 //    } else {
00178 //        cat.setLabel(val.c_str()); 
00179 //    }
00180 //}
00181 
00182 
00183 void CascadeMinimizer::trivialMinimize(const RooAbsReal &nll, RooRealVar &r, int points) const {
00184     double rMin = r.getMin(), rMax = r.getMax(), rStep = (rMax-rMin)/(points-1);
00185     int iMin = -1; double minnll = 0;
00186     for (int i = 0; i < points; ++i) {
00187         double x = rMin + (i+0.5)*rStep;
00188         r.setVal(x);
00189         double y = nll.getVal();
00190         if (iMin == -1 || y < minnll) { minnll = y; iMin = i; }
00191     }
00192     r.setVal( rMin + (iMin+0.5)*rStep );
00193 }
/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc
