Further updates for the minimizers

--- yaml --- svn_rev: 112489 current_ref: refs/heads/CMSSW_7_0_X current_commit: d2a66e3 head_branch: refs/heads/CMSSW_7_0_X migrated_from: v3
cms-sw · Mar 2, 2012 · 3325414 · 3325414
1 parent d66b6c8
commit 3325414
Show file tree

Hide file tree

Showing 8 changed files with 183 additions and 45 deletions.
diff --git a/[refs] b/[refs]
@@ -1,2 +1,2 @@
 ---
-refs/heads/CMSSW_7_0_X: 8d707f24cff1d0e3fdb5bca850b657d54fc2111a
+refs/heads/CMSSW_7_0_X: d2a66e3acfb93dac540352c1a9afb14c7c37895c
diff --git a/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h b/trunk/HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h
@@ -33,6 +33,8 @@ class CascadeMinimizer {
         int          strategy_;
         RooRealVar * poi_; 
 
+        bool improveOnce(int verbose);
+
         /// options configured from command line
         static boost::program_options::options_description options_;
         /// compact information about an algorithm
@@ -51,6 +53,8 @@ class CascadeMinimizer {
         static bool singleNuisFit_;
         /// do first a fit of only the POI
         static bool setZeroPoint_;
+        /// don't do old fallback using robustMinimize 
+        static bool oldFallback_;
 };
 
 #endif
diff --git a/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h b/trunk/HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h
@@ -63,7 +63,7 @@ class SequentialMinimizer {
         bool minimize(double ytol=0, int bigsteps=0, int smallsteps=5);
         bool improve(double ytol=0,  int bigsteps=0, int smallsteps=5);
     private:
-        enum State { Cleared, Ready, Active, Done };
+        enum State { Cleared, Ready, Active, Done, Fixed };
         struct Worker : public OneDimMinimizer {
             Worker() : OneDimMinimizer(), xtol(0), state(Cleared) {}
             Worker(RooAbsReal *nll, RooRealVar *var) : OneDimMinimizer(nll,var), xtol(0), state(Cleared) {}
@@ -72,5 +72,6 @@ class SequentialMinimizer {
         };
         RooAbsReal *nll_;
         std::vector<Worker> workers_;
+        State state_;
 };
 #endif
diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc b/trunk/HiggsAnalysis/CombinedLimit/src/CascadeMinimizer.cc
@@ -13,6 +13,7 @@ bool CascadeMinimizer::preScan_;
 bool CascadeMinimizer::poiOnlyFit_;
 bool CascadeMinimizer::singleNuisFit_;
 bool CascadeMinimizer::setZeroPoint_;
+bool CascadeMinimizer::oldFallback_ = true;
 
 CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi, int initialStrategy) :
     nll_(nll),
@@ -27,33 +28,25 @@ bool CascadeMinimizer::improve(int verbose, bool cascade)
 {
     if (setZeroPoint_) {
         cacheutils::CachingSimNLL *simnll = dynamic_cast<cacheutils::CachingSimNLL *>(&nll_);
-        if (simnll) simnll->setZeroPoint();
+        if (simnll) { 
+            simnll->setZeroPoint();
+        }
     }
     minimizer_.setPrintLevel(verbose-2);  
     minimizer_.setStrategy(strategy_);
-    std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType());
-    std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo());
-    float       nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance());
-    bool outcome = false;
-    if (nominalType == "Sequential") {
-        if (seqmin_.get() == 0) {
-            seqmin_.reset(new SequentialMinimizer(& nll_, mode_ == Unconstrained ? poi_ : 0));
-            seqmin_->minimize(nominalTol);
-        } else {
-            seqmin_->improve(nominalTol);
-        }
-    } else {
-        outcome = nllutils::robustMinimize(nll_, minimizer_, verbose);
-    }
+    bool outcome = improveOnce(verbose-2);
     if (cascade && !outcome && !fallbacks_.empty()) {
+        std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType());
+        std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo());
+        float       nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance());
         if (verbose > 0) std::cerr << "Failed minimization with " << nominalType << "," << nominalAlgo << " and tolerance " << nominalTol << std::endl;
         for (std::vector<Algo>::const_iterator it = fallbacks_.begin(), ed = fallbacks_.end(); it != ed; ++it) {
             ProfileLikelihood::MinimizerSentry minimizerConfig(it->algo, it->tolerance != -1.f ? it->tolerance : nominalTol);
             if (nominalType != ROOT::Math::MinimizerOptions::DefaultMinimizerType() ||
                 nominalAlgo != ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo() ||
                 nominalTol  != ROOT::Math::MinimizerOptions::DefaultTolerance()) {
                 if (verbose > 0) std::cerr << "Will fallback to minimization using " << it->algo << " and tolerance " << it->tolerance << std::endl;
-                outcome = nllutils::robustMinimize(nll_, minimizer_, verbose);
+                outcome = improveOnce(verbose-2);
                 if (outcome) break;
             }
         }
@@ -65,6 +58,28 @@ bool CascadeMinimizer::improve(int verbose, bool cascade)
     return outcome;
 }
 
+bool CascadeMinimizer::improveOnce(int verbose) 
+{
+    std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType());
+    std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo());
+    float       myTol(ROOT::Math::MinimizerOptions::DefaultTolerance());
+    bool outcome = false;
+    if (myType == "Sequential") {
+        if (seqmin_.get() == 0) {
+            seqmin_.reset(new SequentialMinimizer(& nll_, mode_ == Unconstrained ? poi_ : 0));
+            outcome = seqmin_->minimize(myTol);
+        } else {
+            outcome = seqmin_->improve(myTol);
+        }
+    } else if (oldFallback_){
+        outcome = nllutils::robustMinimize(nll_, minimizer_, verbose);
+    } else {
+        int status = minimizer_.minimize(myType.c_str(), myAlgo.c_str());
+        outcome = (status == 0);
+    }
+    return outcome;
+}
+
 bool CascadeMinimizer::minimize(int verbose, bool cascade) 
 {
     minimizer_.setPrintLevel(verbose-2);  
@@ -86,6 +101,7 @@ void CascadeMinimizer::initOptions()
         ("cminSingleNuisFit", "Do first a minimization of each nuisance parameter individually")
         ("cminFallbackAlgo", boost::program_options::value<std::vector<std::string> >(), "Fallback algorithms if the default minimizer fails (can use multiple ones). Syntax is algo[,subalgo][:tolerance]")
         ("cminSetZeroPoint", "Change the reference point of the NLL to be zero during minimization")
+        ("cminOldRobustMinimize", boost::program_options::value<bool>(&oldFallback_)->default_value(oldFallback_), "Use the old 'robustMinimize' logic in addition to the cascade")
         ;
 }
 

diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc b/trunk/HiggsAnalysis/CombinedLimit/src/ProfileLikelihood.cc
@@ -386,7 +386,7 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat
     }
     poi.setConstant(true);
     CascadeMinimizer minim(*nll, CascadeMinimizer::Constrained);
-    if (!minim.improve(verbose-2)) {
+    if (!minim.minimize(verbose-2)) {
         std::cerr << "Initial minimization failed. Aborting." << std::endl;
         return -1;
     } else if (verbose > 0) {
@@ -396,10 +396,14 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat
     std::auto_ptr<RooFitResult> start(minim.save());
     double minnll = nll->getVal(), thisnll = minnll, lastnll = thisnll;
     double rbest = poi.getVal(), rval = rbest;
+    TGraph *points = 0;
     if (verbose) {
         printf("  %-6s  delta(NLL)\n", poi.GetName());
         printf("%8.5f  %8.5f\n", rval, 0.);
         fflush(stdout);
+        points = new TGraph(1); 
+        points->SetName(Form("nll_scan_%g", mass_));
+        points->SetPoint(0, rval, 0);
     }
     while (rval >= tolerance * poi.getMax()) {
         rval *= 0.8;
@@ -412,7 +416,11 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat
             std::cerr << "Minimization failed at " << poi.getVal() <<". exiting the loop" << std::endl;
             return -1;
         } 
-        if (verbose) {  printf("%8.5f  %8.5f\n", rval, thisnll-minnll); fflush(stdout);  }
+        if (verbose) {  
+            printf("%8.5f  %8.5f\n", rval, thisnll-minnll); fflush(stdout);  
+            points->Set(points->GetN()+1);
+            points->SetPoint(points->GetN()-1, rval, thisnll - minnll);
+        }
         if (fabs(lastnll - thisnll) < 7*minimizerToleranceForBF_) {
             std::cout << "This is enough." << std::endl;
             if (thisnll < lastnll) {
@@ -438,6 +446,7 @@ double ProfileLikelihood::significanceBruteForce(RooAbsPdf &pdf, RooAbsData &dat
         }
 #endif
    }
+    if (points) outputFile->WriteTObject(points);
    return (thisnll - minnll);
 }
 

diff --git a/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc b/trunk/HiggsAnalysis/CombinedLimit/src/ProfiledLikelihoodRatioTestStatExt.cc
@@ -14,7 +14,7 @@
 #include "../interface/ProfilingTools.h"
 
 //---- Uncomment this and run with --perfCounters to get statistics of successful and failed fits
-//#define DEBUG_FIT_STATUS
+#define DEBUG_FIT_STATUS
 #ifdef DEBUG_FIT_STATUS
 #define  COUNT_ONE(x) PerfCounter::add(x);
 #else
@@ -147,6 +147,7 @@ ProfiledLikelihoodTestStatOpt::ProfiledLikelihoodTestStatOpt(
 Double_t ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/)
 {
     bool do_debug = runtimedef::get("DEBUG_PLTSO");
+    //static bool do_rescan = runtimedef::get("PLTSO_FULL_RESCAN");
     DBG(DBG_PLTestStat_main, (std::cout << "Being evaluated on " << data.GetName() << std::endl))
 
     // Take snapshot of initial state, to restore it at the end 
@@ -207,6 +208,19 @@ Double_t ProfiledLikelihoodTestStatOpt::Evaluate(RooAbsData& data, RooArgSet& /*
                 thisNLL = nullNLL;
             }
         }
+        /* This piece of debug code was added to see if we were finding a local minimum at zero instead of the global minimum.
+           It doesn't seem to be the case, however  
+        if (do_rescan && fabs(thisNLL - nullNLL) < 0.2 && initialR == 0) {
+            if (do_debug) printf("Doing full rescan. best fit r = %.4f, -lnQ = %.5f\n", bestFitR, thisNLL-nullNLL);
+            for (double rx = 0; rx <= 5; rx += 0.01) {
+                r->setVal(rx); double y = nll_->getVal();
+                if (y < nullNLL) {
+                    printf("%4.2f\t%8.5f\t<<==== ALERT\n", rx, y - nullNLL);
+                } else {
+                    printf("%4.2f\t%8.5f\n", rx, y - nullNLL);
+                }
+            }
+        } */
         if (bestFitR > initialR && oneSided_ == signFlipDef) {
             DBG(DBG_PLTestStat_main, (printf("   fitted signal %7.4f > %7.4f, test statistics will be negative.\n", bestFitR, initialR)))
             std::swap(thisNLL, nullNLL);