Merge pull request #511 from cms-analysis/nckw_freezeUselessParameters

Nckw freeze useless parameters

interface/CascadeMinimizer.h

@@ -36,7 +36,8 @@ class CascadeMinimizer {
         static const boost::program_options::options_description & options() { return options_; }
         void trivialMinimize(const RooAbsReal &nll, RooRealVar &r, int points=100) const ;
         //void collectIrrelevantNuisances(RooAbsCollection &irrelevant) const ;
-        void setAutoBounds(const RooArgSet *pois) ; 
+	bool freezeDiscParams(const bool);
+        void setAutoBounds(const RooArgSet *pois) ;
         void setAutoMax(const RooArgSet *pois) ; 
 	double tolerance() {return defaultMinimizerTolerance_;};
 	std::string algo() {return defaultMinimizerAlgo_;};
@@ -124,6 +125,7 @@ class CascadeMinimizerGlobalConfigs{
 	  RooListProxy allFloatingParameters; 
 	  RooListProxy parametersOfInterest; 
 	  RooListProxy allRooMultiPdfParams;
+	  RooListProxy allRooMultiPdfs;
 
 	  static CascadeMinimizerGlobalConfigs& O(){
 

interface/utils.h

@@ -118,5 +118,6 @@ namespace utils {
     std::vector<std::vector<int> > generateOrthogonalCombinations(const std::vector<int> &vec);
     int countFloating(const RooArgSet &);
     RooArgSet returnAllVars(RooWorkspace *);
+    bool freezeAllDisassociatedRooMultiPdfParameters(RooArgSet multiPdfs, RooArgSet allRooMultiPdfParams, bool freeze=true);
 }
 #endif

src/CascadeMinimizer.cc

@@ -57,6 +57,16 @@ CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi)
 {
 }
 
+bool CascadeMinimizer::freezeDiscParams(const bool freeze)
+{
+    if (runtimedef::get(std::string("MINIMIZER_freezeDiscretes"))) {
+      bool ret =  utils::freezeAllDisassociatedRooMultiPdfParameters((CascadeMinimizerGlobalConfigs::O().allRooMultiPdfs),(CascadeMinimizerGlobalConfigs::O().allRooMultiPdfParams),freeze);
+      return ret;
+    } else {
+      return false;
+    }
+}
+
 void CascadeMinimizer::setAutoBounds(const RooArgSet *pois) 
 {
     poisForAutoBounds_ = pois;
@@ -279,8 +289,22 @@ bool CascadeMinimizer::iterativeMinimize(double &minimumNLL,int verbose, bool ca
    */
 
    // Do A reasonable fit if something changed before 
-   if ( fabs(minimumNLL - nll_.getVal()) > discreteMinTol_ ) improve(verbose,cascade);
+
+   // First freeze all parameters that have nothing to do with the current active pdfs*
+   freezeDiscParams(true);
+
+   //std::cout << " Staring in iterativeMinimize and the minimum NLL so far is  " << minimumNLL << std::endl; 
+   if ( fabs(minimumNLL - nll_.getVal()) > discreteMinTol_ ) { 
+
+     minimizer_.reset(new RooMinimizer(nll_));
+     cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast<cacheutils::CachingSimNLL *>(&nll_) : 0;
+     if (simnll) simnll->setZeroPoint();
+
+     improve(verbose,cascade);
+     //std::cout << " Had to improve further since tolerance is not yet reached   " << nll_.getVal() << std::endl; 
+   }
 
+   // Next remove the POIs and constrained nuisances - this is to set up for the fast loop over the Index combinations
    RooArgSet nuisances = CascadeMinimizerGlobalConfigs::O().allFloatingParameters;
    nuisances.remove(CascadeMinimizerGlobalConfigs::O().allRooMultiPdfParams);
 
@@ -291,7 +315,6 @@ bool CascadeMinimizer::iterativeMinimize(double &minimumNLL,int verbose, bool ca
    if (poi.getSize() >0) frozen.add(poi);
 
    RooStats::RemoveConstantParameters(&frozen);
-
    utils::setAllConstant(frozen,true);
 
    // remake the minimizer   
@@ -315,13 +338,17 @@ bool CascadeMinimizer::iterativeMinimize(double &minimumNLL,int verbose, bool ca
    if (simnll) simnll->clearZeroPoint();
 
    utils::setAllConstant(frozen,false);
-
-   //if (discretesHaveChanged) { 
+
    // Run one last fully floating fit to maintain RooFitResult
    minimizer_.reset(new RooMinimizer(nll_));
    ret = improve(verbose, cascade); 
    //}
    minimumNLL = nll_.getVal();
+   //std::cout << " At the end of iterativeMinimizer, minimum NLL is now " << minimumNLL << std::endl; 
+
+   // unfreeze from *
+   freezeDiscParams(false);
+
    return ret;
 }
 
@@ -340,16 +367,18 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade)
 
     minimizer_->setPrintLevel(verbose-2);  
     minimizer_->setStrategy(strategy_);
-    if (preScan_) minimizer_->minimize("Minuit2","Scan");
-
-
+    //if (preScan_) minimizer_->minimize("Minuit2","Scan");
+
     //if (preFit_ && nuisances != 0) {
 
     RooArgSet nuisances = CascadeMinimizerGlobalConfigs::O().nuisanceParameters;
+
+
     if (preFit_ ) {
         RooArgSet frozen(nuisances);
         RooStats::RemoveConstantParameters(&frozen);
         utils::setAllConstant(frozen,true);
+        freezeDiscParams(true);
 
         minimizer_.reset(new RooMinimizer(nll_));
         minimizer_->setPrintLevel(verbose-2);
@@ -361,8 +390,10 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade)
         minimizer_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
         if (simnll) simnll->clearZeroPoint();
         utils::setAllConstant(frozen,false);
+        freezeDiscParams(false);
         minimizer_.reset(new RooMinimizer(nll_));
     }
+
      // FIXME can be made smarter than this
     /*
     if (mode_ == Unconstrained && poiOnlyFit_) {
@@ -460,6 +491,7 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters,
      Mode 2 -- Full scan over the remaining combinations after mode 1
     */
 
+    //std::cout << " At the start of the looping over the Indeces, minimum NLL is " << minimumNLL << std::endl; 
     // If the barlow-beeston minimisation is being used we can disable it temporarily,
     // saves time if we don't have to call enable/disable on the CMSHistErrorPropagators
     // repeatedly for no purpose
@@ -569,19 +601,33 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters,
 
       if (fitCounter>0) params->assignValueOnly(reallyCleanParameters); // no need to reset from 0'th fit
 
+
+      bool resetNeeded = freezeDiscParams(true);
+
+      // frozen some parameters so need to reset the minimizer
+      if (resetNeeded){
+        minimizer_.reset(new RooMinimizer(nll_));
+        cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast<cacheutils::CachingSimNLL *>(&nll_) : 0;
+        if (simnll) simnll->setZeroPoint();
+      }
+
       // FIXME can be made smarter than this
       if (mode_ == Unconstrained && poiOnlyFit_) {
         trivialMinimize(nll_, *poi_, 200);
       }
 
       ret =  improve(verbose, cascade);
 
+      freezeDiscParams(false);
+
       fitCounter++;
       double thisNllValue = nll_.getVal();
+      //std::cout << " After constrained fit, NLL is " << thisNllValue << std::endl; 
 
       if ( thisNllValue < minimumNLL ){
 		// Now we insert the correction ! 
 	        minimumNLL = thisNllValue;	
+                //std::cout << " .... Found a better fit! hoorah! " << minimumNLL << std::endl; 
     		snap.assignValueOnly(*params);
 		// set the best indeces again
 		for (int id=0;id<numIndeces;id++) {
@@ -626,6 +672,7 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters,
 	((RooCategory*)(pdfCategoryIndeces.at(id)))->setIndex(bestIndeces[id]);	
     } 
     params->assignValueOnly(snap);
+    //std::cout << " After iterations, minimum NLL is " << minimumNLL << std::endl; 
 
     runtimedef::set("MINIMIZER_analytic", currentBarlowBeeston);
     ROOT::Math::MinimizerOptions::SetDefaultStrategy(backupStrategy);
@@ -637,7 +684,7 @@ void CascadeMinimizer::initOptions()
     options_.add_options()
         ("cminPoiOnlyFit",  "Do first a fit floating only the parameter of interest")
         ("cminPreScan",  "Do a scan before first minimization")
-        ("cminPreFit", boost::program_options::value<int>(&preFit_)->default_value(preFit_), "if set to a value N > 0, it will perform a pre-fit with strategy (N-1) with frozen nuisance parameters.")
+        ("cminPreFit", boost::program_options::value<int>(&preFit_)->default_value(preFit_), "if set to a value N > 0, it will perform a pre-fit with strategy (N-1) with frozen constrained nuisance parameters.")
         ("cminApproxPreFitTolerance", boost::program_options::value<double>(&approxPreFitTolerance_)->default_value(approxPreFitTolerance_), "If non-zero, do first a pre-fit with this tolerance (or 10 times the final tolerance, whichever is largest)")
         ("cminApproxPreFitStrategy", boost::program_options::value<int>(&approxPreFitStrategy_)->default_value(approxPreFitStrategy_), "Strategy to use in the pre-fit")
         ("cminSingleNuisFit", "Do first a minimization of each nuisance parameter individually")

src/Combine.cc

@@ -59,6 +59,7 @@
 #include "HiggsAnalysis/CombinedLimit/interface/AsimovUtils.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h"
 #include "HiggsAnalysis/CombinedLimit/interface/ProfilingTools.h"
+#include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h"
 
 #include "HiggsAnalysis/CombinedLimit/interface/Logger.h"
@@ -1228,7 +1229,9 @@ void Combine::addDiscreteNuisances(RooWorkspace *w){
     RooArgSet clients;
     utils::getClients(CascadeMinimizerGlobalConfigs::O().pdfCategories,(w->allPdfs()),clients);
     TIterator *it = clients.createIterator();
-    while (RooAbsArg *arg = (RooAbsArg*)it->Next()) { 
+    clients.Print();
+    while (RooAbsArg *arg = (RooAbsArg*)it->Next()) {
+      (CascadeMinimizerGlobalConfigs::O().allRooMultiPdfs).add(*(dynamic_cast<RooMultiPdf*>(arg)));
       RooAbsPdf *pdf = dynamic_cast<RooAbsPdf*>(arg);
       RooArgSet *pdfPars = pdf->getParameters((const RooArgSet*)0);
       std::auto_ptr<TIterator> iter_v(pdfPars->createIterator());

src/utils.cc

@@ -1,5 +1,6 @@
 #include "HiggsAnalysis/CombinedLimit/interface/utils.h"
 #include "HiggsAnalysis/CombinedLimit/interface/RooSimultaneousOpt.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CascadeMinimizer.h"
 
 #include <cstdio>
 #include <iostream>
@@ -28,6 +29,7 @@
 #include <RooWorkspace.h>
 #include <RooPlot.h>
 #include <RooStats/ModelConfig.h>
+#include <RooStats/RooStatsUtils.h>
 
 #include <boost/algorithm/string/split.hpp>
 #include <boost/algorithm/string/classification.hpp>
@@ -37,6 +39,7 @@
 #include "HiggsAnalysis/CombinedLimit/interface/CloseCoutSentry.h"
 #include "HiggsAnalysis/CombinedLimit/interface/ProfilingTools.h"
 #include "HiggsAnalysis/CombinedLimit/interface/Logger.h"
+#include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
 
 using namespace std;
 
@@ -1007,3 +1010,32 @@ RooArgSet utils::returnAllVars(RooWorkspace *w){
 	args.add(w->allCats());
 	return args;
 }
+
+bool utils::freezeAllDisassociatedRooMultiPdfParameters(RooArgSet multiPdfs, RooArgSet allRooMultiPdfParams, bool freeze){
+
+	RooArgSet *multiPdfParams = (RooArgSet*) allRooMultiPdfParams.Clone();
+
+	// For each multiPdf, get the active pdf and remove its parameters 
+	// from this list of params and then freeze the remaining ones 
+
+        std::auto_ptr<TIterator> iter_p(multiPdfs.createIterator());
+        for (RooAbsArg *P = (RooAbsArg *) iter_p->Next(); P != 0; P = (RooAbsArg *) iter_p->Next()) {
+	  RooMultiPdf *mpdf = dynamic_cast<RooMultiPdf *>(P);
+	  RooAbsPdf *pdf = (RooAbsPdf*)mpdf->getCurrentPdf();
+	  std::cout << " Current active PDF - " << pdf->GetName() <<std::endl;
+	  RooArgSet *pdfPars = pdf->getParameters((const RooArgSet*)0);
+	  RooStats::RemoveConstantParameters(pdfPars); // make sure still to ignore user set constants 
+	  multiPdfParams->remove(*pdfPars);
+
+	} 
+
+	if (multiPdfParams->getSize()>0 ) {
+	 std::cout << " Going to freeze these disacociated params" << std::endl; 
+	 multiPdfParams->Print();
+	 setAllConstant(*multiPdfParams,freeze);
+	 return true;
+	}
+
+	return false;
+}
+