Merge pull request cms-sw#59 from gpetruc/CMSDAS-2013-Kolkata

Welcome to the `development` branch.

BuildFile.xml

@@ -4,6 +4,7 @@
 <use name="roostats"/>
 <use name="histfactory"/>
 <use name="libxml2"/>
+<use name="vdt"/>
 <use   name="boost_program_options"/>
 <use   name="boost_filesystem"/>
 <export>

bin/combine.cpp

@@ -237,8 +237,13 @@ int main(int argc, char **argv) {
 
   if (vm.count("X-fpeMask")) gSystem->SetFPEMask(vm["X-fpeMask"].as<int>());
 
+  // CMSDAS Defaults (you can turn off with --X-rtd <name>=0
+  runtimedef::set("OPTIMIZE_BOUNDS", 1);
+  runtimedef::set("ADDNLL_RECURSIVE", 1);
+  runtimedef::set("ADDNLL_GAUSSNLL", 1);
+  runtimedef::set("ADDNLL_HISTNLL", 1);
+  runtimedef::set("TMCSO_AdaptivePseudoAsimov", 1);
 
-  // if you have libraries, it's time to load them now
   for (vector<string>::const_iterator rtdp = runtimeDefines.begin(), endrtdp = runtimeDefines.end(); rtdp != endrtdp; ++rtdp) {
     std::string::size_type idx = rtdp->find('=');
     if (idx == std::string::npos) {

interface/CachingNLL.h

@@ -13,6 +13,11 @@
 #include <RooRealVar.h>
 #include <RooSimultaneous.h>
 #include <RooGaussian.h>
+#include <RooProduct.h>
+#include "../interface/SimpleGaussianConstraint.h"
+#include <boost/ptr_container/ptr_vector.hpp>
+
+class RooMultiPdf;
 
 // Part zero: ArgSet checker
 namespace cacheutils {
@@ -55,20 +60,38 @@ class CachingPdf {
     public:
         CachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ;
         CachingPdf(const CachingPdf &other) ;
-        ~CachingPdf() ;
+        virtual ~CachingPdf() ;
         const std::vector<Double_t> & eval(const RooAbsData &data) ;
         const RooAbsReal *pdf() const { return pdf_; }
         void  setDataDirty() { lastData_ = 0; }
-    private:
+    protected:
         const RooArgSet *obs_;
         RooAbsReal *pdfOriginal_;
         RooArgSet  pdfPieces_;
         RooAbsReal *pdf_;
         const RooAbsData *lastData_;
         ValuesCache cache_;
-        void realFill_(const RooAbsData &data, std::vector<Double_t> &values) ;
+        std::vector<uint8_t> nonZeroW_;
+        unsigned int         nonZeroWEntries_;
+        virtual void newData_(const RooAbsData &data) ;
+        virtual void realFill_(const RooAbsData &data, std::vector<Double_t> &values) ;
 };
 
+template <typename PdfT, typename VPdfT> 
+class OptimizedCachingPdfT : public CachingPdf {
+    public:
+        OptimizedCachingPdfT(RooAbsReal *pdf, const RooArgSet *obs) :
+            CachingPdf(pdf,obs), vpdf_(0) {}
+        OptimizedCachingPdfT(const OptimizedCachingPdfT &other) : 
+            CachingPdf(other), vpdf_(0) {}
+        virtual ~OptimizedCachingPdfT() { delete vpdf_; }
+    protected:
+        virtual void realFill_(const RooAbsData &data, std::vector<Double_t> &values) ;
+        virtual void newData_(const RooAbsData &data) ;
+        VPdfT *vpdf_;
+};
+
+
 class CachingAddNLL : public RooAbsReal {
     public:
         CachingAddNLL(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data) ;
@@ -88,16 +111,20 @@ class CachingAddNLL : public RooAbsReal {
         RooSetProxy & params() { return params_; }
     private:
         void setup_();
+        void addPdfs_(RooAddPdf *addpdf, bool recursive, const RooArgList & basecoeffs) ;
         RooAbsPdf *pdf_;
         RooSetProxy params_;
         const RooAbsData *data_;
         std::vector<Double_t>  weights_;
         double               sumWeights_;
         mutable std::vector<RooAbsReal*> coeffs_;
-        mutable std::vector<CachingPdf>  pdfs_;
+        mutable boost::ptr_vector<CachingPdf>  pdfs_;
+        mutable boost::ptr_vector<RooAbsReal>  prods_;
         mutable std::vector<RooAbsReal*> integrals_;
+        mutable std::vector<std::pair<const RooMultiPdf*,CachingPdf*> > multiPdfs_;
         mutable std::vector<Double_t> partialSum_;
-        mutable bool isRooRealSum_;
+        mutable std::vector<Double_t> workingArea_;
+        mutable bool isRooRealSum_, fastExit_;
         double zeroPoint_;
 };
 
@@ -117,18 +144,9 @@ class CachingSimNLL  : public RooAbsReal {
         static void setNoDeepLogEvalError(bool noDeep) { noDeepLEE_ = noDeep; }
         void setZeroPoint() ; 
         void clearZeroPoint() ;
+        static void forceUnoptimizedConstraints() { optimizeContraints_ = false; }
         friend class CachingAddNLL;
     private:
-        class SimpleGaussianConstraint : public RooGaussian {
-            public:
-                SimpleGaussianConstraint(const RooGaussian &g) : RooGaussian(g, "") {}
-                double getLogValFast() const { 
-                    Double_t arg = x - mean;  
-                    Double_t sig = sigma ;
-                    return -0.5*arg*arg/(sig*sig);
-                }
-        };
-
         void setup_();
         RooSimultaneous   *pdfOriginal_;
         const RooAbsData  *dataOriginal_;
@@ -138,11 +156,13 @@ class CachingSimNLL  : public RooAbsReal {
         std::auto_ptr<RooSimultaneous>  factorizedPdf_;
         std::vector<RooAbsPdf *>        constrainPdfs_;
         std::vector<SimpleGaussianConstraint *>  constrainPdfsFast_;
+        std::vector<bool>                        constrainPdfsFastOwned_;
         std::vector<CachingAddNLL*>     pdfs_;
         std::auto_ptr<TList>            dataSets_;
         std::vector<RooDataSet *>       datasets_;
         static bool noDeepLEE_;
         static bool hasError_;
+        static bool optimizeContraints_;
         std::vector<double> constrainZeroPoints_;
         std::vector<double> constrainZeroPointsFast_;
 };

interface/Combine.h

@@ -71,6 +71,7 @@ class Combine {
   bool makeTempDir_;
   bool rebuildSimPdf_;
   bool optSimPdf_;
+  bool noMCbonly_;
 
   static TTree *tree_;
 };

interface/FastTemplate.h

@@ -4,31 +4,36 @@
 #include <TH1.h>
 #include <TH2.h>
 #include <algorithm>
+#include <vector>
 
 class FastTemplate {
     public:
         typedef double T;
-        FastTemplate() : size_(0), values_(0) {}
-        FastTemplate(unsigned int size) : size_(size), values_(new T[size_]) {}
-        FastTemplate(const FastTemplate &other) : size_(other.size_), values_(size_ ? new T[size_] :  0) { if (size_) CopyValues(other); }
-        FastTemplate(const TH1 &other) : size_(other.GetNbinsX()), values_(new T[size_]) { CopyValues(other); }
-        FastTemplate(const TH2 &other) : size_(other.GetNbinsX()*other.GetNbinsY()), values_(new T[size_]) { CopyValues(other); }
+        typedef std::vector<T> AT;
+        FastTemplate() : size_(0), values_() {}
+        FastTemplate(unsigned int size) : size_(size), values_(size_) {}
+        FastTemplate(const FastTemplate &other) : size_(other.size()), values_(other.values_) {}
+        FastTemplate(const TH1 &other) : size_(other.GetNbinsX()), values_(size_) { CopyValues(other); }
+        FastTemplate(const TH2 &other) : size_(other.GetNbinsX()*other.GetNbinsY()), values_(size_) { CopyValues(other); }
         FastTemplate & operator=(const FastTemplate &other) { 
-            if (size_ != other.size_) { 
-                delete [] values_; size_ = other.size_; values_ = new T[size_];
+            if (&other != this) {
+                size_ = other.size_;
+                values_ = other.values_;
             }
-            CopyValues(other); return *this; 
+            return *this; 
         }
         FastTemplate & operator=(const TH1 &other) { 
-            if (size_ != unsigned(other.GetNbinsX())) { 
-                delete [] values_; size_ = other.GetNbinsX(); values_ = new T[size_];
+            if (size() != unsigned(other.GetNbinsX())) { 
+                size_ = other.GetNbinsX();
+                values_.resize(size_);
             }
             CopyValues(other); return *this;  
         }
-        ~FastTemplate() { delete [] values_; }
+        ~FastTemplate() {  }
         void Resize(unsigned int newsize) {
-            if (newsize != size_) {
-                delete [] values_; size_ = newsize; values_ = new T[size_];
+            if (newsize != size()) {
+                size_ = newsize;
+                values_.resize(size_);
             }
         }
         T Integral() const ;
@@ -39,6 +44,10 @@ class FastTemplate {
         void CopyValues(const TH2 &other) ;
         T & operator[](unsigned int i) { return values_[i]; }
         const T & operator[](unsigned int i) const { return values_[i]; }
+        /// return the full size of the template
+        const unsigned int fullsize() const { return values_.size(); }
+        /// return the active size of the template (can be less than the full size if the SetActiveSize
+        /// has been used to inform the code that only the first N bins are not empty)
         const unsigned int size() const { return size_; }
 
         /// *this = log(*this) 
@@ -54,40 +63,48 @@ class FastTemplate {
         /// Does this += x * (diff + (sum)*y)
         void Meld(const FastTemplate & diff, const FastTemplate & sum, T x, T y) ;
         /// protect from underflows (*this = max(*this, minimum));
-        void CropUnderflows(T minimum=1e-9);
+        void CropUnderflows(T minimum=1e-9, bool activebinsonly=true);
+
+        /// Tell the code that only the first N bins of the template are non-empty,
+        /// and so that only those have to be considered when doing operations
+        void SetActiveSize(unsigned int size) { size_ = size; }
 
         void Dump() const ;
     protected:
-        unsigned int size_; 
-        T *values_;
+        unsigned int size_;
+        AT values_;
 };
 class FastHisto : public FastTemplate {
     public:
-        FastHisto() : FastTemplate(), binEdges_(0), binWidths_(0) {}
+        FastHisto() : FastTemplate(), binEdges_(), binWidths_() {}
         FastHisto(const TH1 &hist) ;
         FastHisto(const FastHisto &other) ;
         FastHisto & operator=(const FastHisto &other) { 
-            if (size_ != other.size_) {
-                FastHisto fh(other);
-                swap(fh);
+            if (size() != other.size()) {
+                size_ = other.size_;
+                values_    = other.values_;
+                binWidths_ = other.binWidths_;
+                binEdges_  = other.binEdges_;
             } else CopyValues(other); 
             return *this; 
         }
         FastHisto & operator=(const TH1 &other) { 
-            if (size_ != unsigned(other.GetNbinsX())) { 
+            if (size() != unsigned(other.GetNbinsX())) { 
                 FastHisto fh(other);
                 swap(fh);
             } else CopyValues(other); 
-            CopyValues(other); return *this;  
+            return *this;  
         }
-        ~FastHisto() { delete [] binEdges_; delete [] binWidths_; }
+        ~FastHisto() {  }
         void swap(FastHisto &other) {
             std::swap(size_, other.size_);
             std::swap(values_, other.values_);
             std::swap(binWidths_, other.binWidths_);
             std::swap(binEdges_, other.binEdges_);
         }
         T GetAt(const T &x) const ;
+        int FindBin(const T &x) const ;
+        const T & GetBinContent(int bin) const { return values_[bin]; }
         T IntegralWidth() const ;
         void Normalize() {
             T sum = IntegralWidth();
@@ -96,34 +113,41 @@ class FastHisto : public FastTemplate {
 
         void Dump() const ;
     private:
-        T *binEdges_;
-        T *binWidths_;
+        AT binEdges_;
+        AT binWidths_;
 
 };
 class FastHisto2D : public FastTemplate {
     public:
-        FastHisto2D() : FastTemplate(), binX_(0), binY_(0), binEdgesX_(0), binEdgesY_(0), binWidths_(0) {}
+        FastHisto2D() : FastTemplate(), binX_(0), binY_(0), binEdgesX_(), binEdgesY_(), binWidths_() {}
         FastHisto2D(const TH2 &hist, bool normXonly=false) ;
         FastHisto2D(const FastHisto2D &other) ;
         FastHisto2D & operator=(const FastHisto2D &other) { 
-            if (binX_ != other.binY_ || binY_ != other.binY_) {
-                FastHisto2D fh(other);
-                swap(fh);
+            if (binX() != other.binX() || binY() != other.binY()) {
+                size_      = other.size_;
+                values_    = other.values_;
+                binWidths_ = other.binWidths_;
+                binEdgesX_ = other.binEdgesX_;
+                binEdgesY_ = other.binEdgesY_;
+                binX_      = other.binX_;
+                binY_      = other.binY_;
             } else CopyValues(other); 
             return *this; 
         }
-        ~FastHisto2D() { delete [] binEdgesX_; delete [] binEdgesY_; delete [] binWidths_; }
+        ~FastHisto2D() {  }
         void swap(FastHisto2D &other) {
+            std::swap(size_, other.size_);
             std::swap(binX_, other.binX_);
             std::swap(binY_, other.binY_);
-            std::swap(size_, other.size_);
             std::swap(values_, other.values_);
             std::swap(binWidths_, other.binWidths_);
             std::swap(binEdgesX_, other.binEdgesX_);
             std::swap(binEdgesY_, other.binEdgesY_);
         }
         T GetAt(const T &x, const T &y) const ;
         T IntegralWidth() const ;
+        unsigned int binX() const { return binX_; }
+        unsigned int binY() const { return binY_; }
         void Normalize() {
             T sum = IntegralWidth();
             if (sum > 0) Scale(1.0f/sum);
@@ -134,9 +158,9 @@ class FastHisto2D : public FastTemplate {
         void Dump() const ;
     private:
         unsigned int binX_, binY_;
-        T *binEdgesX_;
-        T *binEdgesY_;
-        T *binWidths_;
+        AT binEdgesX_;
+        AT binEdgesY_;
+        AT binWidths_;
 
 };
 

interface/FitterAlgoBase.h

@@ -56,6 +56,9 @@ class FitterAlgoBase : public LimitAlgo {
   RooFitResult *doFit(RooAbsPdf &pdf, RooAbsData &data, const RooArgList &rs, const RooCmdArg &constrain, bool doHesse=true, int ndim=1,bool reuseNLL=false, bool saveFitResult=true) ;
   double findCrossing(CascadeMinimizer &minim, RooAbsReal &nll, RooRealVar &r, double level, double rStart, double rBound) ;
   double findCrossingNew(CascadeMinimizer &minim, RooAbsReal &nll, RooRealVar &r, double level, double rStart, double rBound) ;
+
+  void optimizeBounds(const RooWorkspace *w, const RooStats::ModelConfig *mc) ;
+  void restoreBounds(const RooWorkspace *w, const RooStats::ModelConfig *mc) ;
 };
 
 

interface/RooMinimizerOpt.h

@@ -13,15 +13,42 @@ class RooMinimizerOpt : public RooMinimizer {
     public:
         RooMinimizerOpt(RooAbsReal& function) ;
         Double_t edm();
+        Int_t minimize(const char* type, const char* alg=0) ;
+        Int_t improve() ;
+        Int_t migrad() ;
+        Int_t hesse() ;
+        Int_t minos() ;
+        Int_t minos(const RooArgSet& minosParamList) ;
+
 };
 
 class RooMinimizerFcnOpt : public RooMinimizerFcn {
     public: 
         RooMinimizerFcnOpt(RooAbsReal *funct, RooMinimizer *context,  bool verbose = false);
+        RooMinimizerFcnOpt(const RooMinimizerFcnOpt &other) ;
         virtual ROOT::Math::IBaseFunctionMultiDim* Clone() const;
+        Bool_t Synchronize(std::vector<ROOT::Fit::ParameterSettings>& parameters, Bool_t optConst, Bool_t verbose);
+        void initStdVects() const ;
     protected:
         virtual double DoEval(const double * x) const;
         mutable std::vector<RooRealVar *> _vars;
+        mutable std::vector<double>       _vals;
+        mutable std::vector<bool  >       _hasOptimzedBounds;
+        struct OptBound {
+            double softMin, softMax, hardMin, hardMax;
+            double transform(double x) const {
+                if (x < softMin) {
+                    double dx = softMin-x; // > 0
+                    return hardMin + ( (softMin-hardMin) + dx ) * std::exp ( -2*dx/(softMin-hardMin) );
+                } else if (x > softMax) { 
+                    double dx = x-softMax;
+                    return hardMax - ( (hardMax-softMax) + dx ) * std::exp ( -2*dx/(hardMax-softMax) );
+                } else {
+                    return x;
+                }
+            }
+        };
+        mutable std::vector<OptBound> _optimzedBounds;
 };
 
 #endif