Do analytical Barlow Beeston before likelihood evaluation

Doing the analytical Barlow Beeston before evaluating the likelihood has
two big advantages, instead of doing it while evaluating the
CMSHistSum/CMSHistErrorPropagator:

1. We don't have to hack the dirty state propagation. This was
   necessary before, because it is not allowed to change a RooRealVar
   server during the evaluation of the client. Also, the hack relies on
   the implementation details of the CachingSimNLL (evaluating main pdf
   before constraints), so it's also a problem for using RooFit regular
   NLL code path.

2. If the analytical minimization and the likelihood evaluation
   factorize, we don't need to worry how to differentiate through the
   analytical minimization and the likelihood evaluation together.

Author: guitargeek

interface/CMSHistErrorPropagator.h

@@ -30,7 +30,6 @@ class CMSHistErrorPropagator : public RooAbsReal {
     std::vector<double> x2;
     std::vector<double> res;
     std::vector<double> gobs;
-    std::set<RooAbsArg*> dirty_prop;
     std::vector<RooRealVar*> push_res;
   };
 public:
@@ -77,7 +76,9 @@ class CMSHistErrorPropagator : public RooAbsReal {
 
   friend class CMSHistV<CMSHistErrorPropagator>;
 
- protected:
+  void runBarlowBeeston() const;
+
+protected:
   RooRealProxy x_;
   RooListProxy funcs_;
   RooListProxy coeffs_;
@@ -109,8 +110,6 @@ class CMSHistErrorPropagator : public RooAbsReal {
   void initialize() const;
   void updateCache(int eval = 1) const;
 
-  void runBarlowBeeston() const;
-
 
  private:
   ClassDefOverride(CMSHistErrorPropagator,1)

interface/CMSHistSum.h

@@ -30,7 +30,6 @@ class CMSHistSum : public RooAbsReal {
     std::vector<double> x2;
     std::vector<double> res;
     std::vector<double> gobs;
-    std::set<RooAbsArg*> dirty_prop;
     std::vector<RooRealVar*> push_res;
   };
 public:
@@ -80,7 +79,9 @@ class CMSHistSum : public RooAbsReal {
 
   void injectExternalMorph(int idx, CMSExternalMorph& morph);
 
- protected:
+  void runBarlowBeeston() const;
+
+protected:
   RooRealProxy x_;
 
   RooListProxy morphpars_;
@@ -143,9 +144,6 @@ class CMSHistSum : public RooAbsReal {
   void updateCache() const;
   inline double smoothStepFunc(double x, int const& ip) const;
 
-
-  void runBarlowBeeston() const;
-
   void updateMorphs() const;
 
 

interface/CachingNLL.h

@@ -129,6 +129,7 @@ class CachingAddNLL : public RooAbsReal {
         void updateZeroPoint() { clearZeroPoint(); setZeroPoint(); }
         void propagateData();
         void setAnalyticBarlowBeeston(bool flag);
+        void runAnalyticBarlowBeeston();
         /// note: setIncludeZeroWeights(true) won't have effect unless you also re-call setData
         virtual void  setIncludeZeroWeights(bool includeZeroWeights) ;
         RooSetProxy & params() { return params_; }

src/CMSHistErrorPropagator.cc

@@ -99,7 +99,6 @@ void CMSHistErrorPropagator::initialize() const {
   initialized_ = true;
 }
 
-
 void CMSHistErrorPropagator::updateCache(int eval) const {
   initialize();
 
@@ -153,7 +152,6 @@ void CMSHistErrorPropagator::updateCache(int eval) const {
 
 
   if (!binsentry_.good() || eval != last_eval_) {
-    runBarlowBeeston();
     // bintypes might have size == 0 if we never ran setupBinPars()
     for (unsigned j = 0; j < bintypes_.size(); ++j) {
       cache_[j] = valsum_[j];
@@ -192,9 +190,7 @@ void CMSHistErrorPropagator::updateCache(int eval) const {
 }
 
 void CMSHistErrorPropagator::runBarlowBeeston() const {
-  if (!bb_.init) return;
-  RooAbsArg::setDirtyInhibit(true);
-
+  updateCache();
   const unsigned n = bb_.use.size();
   for (unsigned j = 0; j < n; ++j) {
     bb_.dat[j] = data_[bb_.use[j]];
@@ -215,10 +211,6 @@ void CMSHistErrorPropagator::runBarlowBeeston() const {
   for (unsigned j = 0; j < n; ++j) {
     if (toterr_[bb_.use[j]] > 0.) bb_.push_res[j]->setVal(bb_.res[j]);
   }
-  RooAbsArg::setDirtyInhibit(false);
-  for (RooAbsArg *arg : bb_.dirty_prop) {
-    arg->setValueDirty();
-  }
 }
 
 void CMSHistErrorPropagator::setAnalyticBarlowBeeston(bool flag) const {
@@ -240,7 +232,6 @@ void CMSHistErrorPropagator::setAnalyticBarlowBeeston(bool flag) const {
     bb_.x2.clear();
     bb_.res.clear();
     bb_.gobs.clear();
-    bb_.dirty_prop.clear();
     bb_.push_res.clear();
     bb_.init = false;
   }
@@ -254,7 +245,6 @@ void CMSHistErrorPropagator::setAnalyticBarlowBeeston(bool flag) const {
             // std::cout << "Skipping " << this << " " << this->GetName() << "\n";
           } else {
             // std::cout << "Adding " << arg << " " << arg->GetName() << "\n";
-            bb_.dirty_prop.insert(arg);
             auto as_gauss = dynamic_cast<RooGaussian*>(arg);
             if (as_gauss) {
               auto gobs = dynamic_cast<RooAbsReal*>(as_gauss->findServer(TString(vbinpars_[j][0]->GetName())+"_In"));

src/CMSHistSum.cc

@@ -336,7 +336,6 @@ inline double CMSHistSum::smoothStepFunc(double x, int const& ip) const {
   return 0.125 * xnorm * (xnorm2 * (3. * xnorm2 - 10.) + 15);
 }
 
-
 void CMSHistSum::updateCache() const {
   initialize();
 
@@ -379,10 +378,6 @@ void CMSHistSum::updateCache() const {
 
 
   if (!binsentry_.good()) {
-    #if HFVERBOSE > 0
-      std::cout << "Calling runBarlowBeeston\n";
-    #endif
-    runBarlowBeeston();
     // bintypes might have size == 0 if we never ran setupBinPars()
     #if HFVERBOSE > 0
       std::cout << "Assigning bin shifts\n";
@@ -418,9 +413,7 @@ void CMSHistSum::updateCache() const {
 }
 
 void CMSHistSum::runBarlowBeeston() const {
-  if (!bb_.init) return;
-  RooAbsArg::setDirtyInhibit(true);
-
+  updateCache();
   const unsigned n = bb_.use.size();
   for (unsigned j = 0; j < n; ++j) {
     bb_.dat[j] = data_[bb_.use[j]];
@@ -441,10 +434,6 @@ void CMSHistSum::runBarlowBeeston() const {
   for (unsigned j = 0; j < n; ++j) {
     if (toterr_[bb_.use[j]] > 0.) bb_.push_res[j]->setVal(bb_.res[j]);
   }
-  RooAbsArg::setDirtyInhibit(false);
-  for (RooAbsArg *arg : bb_.dirty_prop) {
-    arg->setValueDirty();
-  }
 }
 
 void CMSHistSum::setAnalyticBarlowBeeston(bool flag) const {
@@ -466,7 +455,6 @@ void CMSHistSum::setAnalyticBarlowBeeston(bool flag) const {
     bb_.x2.clear();
     bb_.res.clear();
     bb_.gobs.clear();
-    bb_.dirty_prop.clear();
     bb_.push_res.clear();
     bb_.init = false;
   }
@@ -480,7 +468,6 @@ void CMSHistSum::setAnalyticBarlowBeeston(bool flag) const {
             // std::cout << "Skipping " << this << " " << this->GetName() << "\n";
           } else {
             // std::cout << "Adding " << arg << " " << arg->GetName() << "\n";
-            bb_.dirty_prop.insert(arg);
             auto as_gauss = dynamic_cast<RooGaussian*>(arg);
             if (as_gauss) {
               auto gobs = dynamic_cast<RooAbsReal*>(as_gauss->findServer(TString(vbinpars_[j][0]->GetName())+"_In"));

src/CachingNLL.cc

@@ -595,6 +595,9 @@ cacheutils::CachingAddNLL::evaluate() const
 #ifdef DEBUG_CACHE
     PerfCounter::add("CachingAddNLL::evaluate called");
 #endif
+    // The very first thing we do before any evaluation: run the analytical
+    // minimization of Barlow-Beeston nuisance parameters.
+    const_cast<CachingAddNLL&>(*this).runAnalyticBarlowBeeston();
 
     std::fill( partialSum_.begin(), partialSum_.end(), 0.0 );
 
@@ -809,6 +812,15 @@ void cacheutils::CachingAddNLL::setAnalyticBarlowBeeston(bool flag) {
   }
 }
 
+void cacheutils::CachingAddNLL::runAnalyticBarlowBeeston() {
+  for (auto* hist : histErrorPropagators_) {
+    hist->runBarlowBeeston();
+  }
+  for (auto* hist : histSums_) {
+    hist->runBarlowBeeston();
+  }
+}
+
 cacheutils::CachingSimNLL::CachingSimNLL(RooSimultaneous *pdf, RooAbsData *data, const RooArgSet *nuis) :
     pdfOriginal_(pdf),
     dataOriginal_(data),
@@ -1007,6 +1019,17 @@ cacheutils::CachingSimNLL::evaluate() const
 #ifdef DEBUG_CACHE
     PerfCounter::add("CachingSimNLL::evaluate called");
 #endif
+
+    // The very first thing we do before any evaluation: run the analytical
+    // minimization of Barlow-Beeston nuisance parameters.
+    for (size_t i = 0; i < pdfs_.size(); ++i) {
+      if (!pdfs_[i])
+        continue;
+      if (!channelMasks_.empty() && channelMasks_[i]->getVal() != 0.)
+        continue;
+      pdfs_[i]->runAnalyticBarlowBeeston();
+    }
+
     static bool gentleNegativePenalty_ = runtimedef::get("GENTLE_LEE");
     DefaultAccumulator<double> ret = 0;
     for (std::size_t idx = 0; idx < pdfs_.size(); ++idx) {