From 47f6feab7ad8f553cfd5f42eb7e122d49bef6446 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 15 Oct 2019 15:37:35 +0200 Subject: [PATCH 01/16] construction of the semianalytic minimizer --- interface/RooFitResult.h | 207 ++++++ interface/RooMinimizerFcnSemiAnalytic.h | 118 +++ interface/RooMinimizerSemiAnalytic.h | 140 ++++ src/RooMinimizerFcnSemiAnalytic.cxx | 716 ++++++++++++++++++ src/RooMinimizerSemiAnalytic.cc | 942 ++++++++++++++++++++++++ src/classes.h | 3 + src/classes_def.xml | 2 + 7 files changed, 2128 insertions(+) create mode 100644 interface/RooFitResult.h create mode 100644 interface/RooMinimizerFcnSemiAnalytic.h create mode 100644 interface/RooMinimizerSemiAnalytic.h create mode 100644 src/RooMinimizerFcnSemiAnalytic.cxx create mode 100644 src/RooMinimizerSemiAnalytic.cc diff --git a/interface/RooFitResult.h b/interface/RooFitResult.h new file mode 100644 index 00000000000..47224da1ec0 --- /dev/null +++ b/interface/RooFitResult.h @@ -0,0 +1,207 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * File: $Id: RooFitResult.h,v 1.28 2007/05/11 09:11:30 verkerke Exp $ + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2005, Regents of the University of California * + * and Stanford University. All rights reserved. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ +#ifndef ROO_FIT_RESULT +#define ROO_FIT_RESULT + +#include "RooAbsArg.h" +#include "RooPrintable.h" +#include "RooDirItem.h" +#include "RooArgList.h" + +#include "RVersion.h" +#include "TMatrixFfwd.h" +#include "TMatrixDSym.h" +#include "TRootIOCtor.h" + +#include +#include +#include + +class RooArgSet ; +class RooAbsPdf ; +class RooPlot; +class TObject ; +class TH2 ; +typedef RooArgSet* pRooArgSet ; + +class RooFitResult : public TNamed, public RooPrintable, public RooDirItem { +public: + + // Constructors, assignment etc. + RooFitResult(const char* name=0, const char* title=0) ; + RooFitResult(const RooFitResult& other) ; + virtual TObject* Clone(const char* newname = 0) const { + RooFitResult* r = new RooFitResult(*this) ; + if (newname && *newname) r->SetName(newname) ; + return r ; + } + virtual TObject* clone() const { return new RooFitResult(*this); } + virtual ~RooFitResult() ; + + static RooFitResult* lastMinuitFit(const RooArgList& varList=RooArgList()) ; + + static RooFitResult *prefitResult(const RooArgList ¶mList); + + // Printing interface (human readable) + virtual void printValue(std::ostream& os) const ; + virtual void printName(std::ostream& os) const ; + virtual void printTitle(std::ostream& os) const ; + virtual void printClassName(std::ostream& os) const ; + virtual void printArgs(std::ostream& os) const ; + void printMultiline(std::ostream& os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const ; + + inline virtual void Print(Option_t *options= 0) const { + // Printing interface + printStream(defaultPrintStream(),defaultPrintContents(options),defaultPrintStyle(options)); + } + + virtual Int_t defaultPrintContents(Option_t* opt) const ; + virtual StyleOption defaultPrintStyle(Option_t* opt) const ; + + RooAbsPdf* createHessePdf(const RooArgSet& params) const ; + + // Accessors + inline Int_t status() const { + // Return MINUIT status code + return _status ; + } + + inline UInt_t numStatusHistory() const { return _statusHistory.size() ; } + Int_t statusCodeHistory(UInt_t icycle) const ; + const char* statusLabelHistory(UInt_t icycle) const ; + + inline Int_t covQual() const { + // Return MINUIT quality code of covariance matrix + return _covQual ; + } + inline Int_t numInvalidNLL() const { + // Return number of NLL evaluations with problems + return _numBadNLL ; + } + inline Double_t edm() const { + // Return estimated distance to minimum + return _edm ; + } + inline Double_t minNll() const { + // Return minimized -log(L) value + return _minNLL ; + } + inline const RooArgList& constPars() const { + // Return list of constant parameters + return *_constPars ; + } + inline const RooArgList& floatParsInit() const { + // Return list of floating parameters before fit + return *_initPars ; + } + inline const RooArgList& floatParsFinal() const { + // Return list of floarting parameters after fit + return *_finalPars ; + } + + TH2* correlationHist(const char* name = "correlation_matrix") const ; + + Double_t correlation(const RooAbsArg& par1, const RooAbsArg& par2) const { + // Return correlation between par1 and par2 + return correlation(par1.GetName(),par2.GetName()) ; + } + const RooArgList* correlation(const RooAbsArg& par) const { + // Return pointer to list of correlations of all parameters with par + return correlation(par.GetName()) ; + } + + Double_t correlation(const char* parname1, const char* parname2) const ; + const RooArgList* correlation(const char* parname) const ; + + + const TMatrixDSym& covarianceMatrix() const ; + const TMatrixDSym& correlationMatrix() const ; + TMatrixDSym reducedCovarianceMatrix(const RooArgList& params) const ; + TMatrixDSym conditionalCovarianceMatrix(const RooArgList& params) const ; + + + // Global correlation accessors + Double_t globalCorr(const RooAbsArg& par) { return globalCorr(par.GetName()) ; } + Double_t globalCorr(const char* parname) ; + const RooArgList* globalCorr() ; + + + // Add objects to a 2D plot + inline RooPlot *plotOn(RooPlot *frame, const RooAbsArg &par1, const RooAbsArg &par2, + const char *options= "ME") const { + // Plot error ellipse in par1 and par2 on frame + return plotOn(frame,par1.GetName(),par2.GetName(),options); + } + RooPlot *plotOn(RooPlot *plot, const char *parName1, const char *parName2, + const char *options= "ME") const; + + // Generate random perturbations of the final parameters using the covariance matrix + const RooArgList& randomizePars() const; + + Bool_t isIdentical(const RooFitResult& other, Double_t tol=5e-5, Double_t tolCorr=1e-4, Bool_t verbose=kTRUE) const ; + + void SetName(const char *name) ; + void SetNameTitle(const char *name, const char* title) ; + +protected: + + friend class RooMinuit ; + friend class RooMinimizer ; + friend class RooMinimizerSemiAnalytic; // :( why with friendship? + + void setCovarianceMatrix(TMatrixDSym& V) ; + void setConstParList(const RooArgList& list) ; + void setInitParList(const RooArgList& list) ; + void setFinalParList(const RooArgList& list) ; + inline void setMinNLL(Double_t val) { _minNLL = val ; } + inline void setEDM(Double_t val) { _edm = val ; } + inline void setStatus(Int_t val) { _status = val ; } + inline void setCovQual(Int_t val) { _covQual = val ; } + inline void setNumInvalidNLL(Int_t val) { _numBadNLL=val ; } + void fillCorrMatrix() ; + void fillCorrMatrix(const std::vector& globalCC, const TMatrixDSym& corrs, const TMatrixDSym& covs) ; + void fillLegacyCorrMatrix() const ; + void fillPrefitCorrMatrix(); + void setStatusHistory(std::vector >& hist) { _statusHistory = hist ; } + + Double_t correlation(Int_t row, Int_t col) const; + Double_t covariance(Int_t row, Int_t col) const; + + Int_t _status ; // MINUIT status code + Int_t _covQual ; // MINUIT quality code of covariance matrix + Int_t _numBadNLL ; // Number calls with bad (zero,negative) likelihood + Double_t _minNLL ; // NLL at minimum + Double_t _edm ; // Estimated distance to minimum + RooArgList* _constPars ; // List of constant parameters + RooArgList* _initPars ; // List of floating parameters with initial values + RooArgList* _finalPars ; // List of floating parameters with final values + + mutable RooArgList* _globalCorr ; //! List of global correlation coefficients + mutable TList _corrMatrix ; //! Correlation matrix (list of RooArgLists) + + mutable RooArgList *_randomPars; //! List of floating parameters with most recent random perturbation applied + mutable TMatrixF* _Lt; //! triangular matrix used for generate random perturbations + + TMatrixDSym* _CM ; // Correlation matrix + TMatrixDSym* _VM ; // Covariance matrix + TVectorD* _GC ; // Global correlation coefficients + + std::vector > _statusHistory ; // History of status codes + + ClassDef(RooFitResult,5) // Container class for fit result +}; + +#endif diff --git a/interface/RooMinimizerFcnSemiAnalytic.h b/interface/RooMinimizerFcnSemiAnalytic.h new file mode 100644 index 00000000000..a885fd27cfc --- /dev/null +++ b/interface/RooMinimizerFcnSemiAnalytic.h @@ -0,0 +1,118 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifndef ROO_MINIMIZER_FCN_SEMIANALYTIC +#define ROO_MINIMIZER_FCN_SEMIANALYTIC + +#include "Math/IFunction.h" +#include "Math/MultiNumGradFunction.h" +#include "Fit/ParameterSettings.h" +#include "Fit/FitResult.h" + +#include "TMatrixDSym.h" + +#include "RooAbsReal.h" +#include "RooArgList.h" + +#include +#include +#include + +class RooMinimizerSemiAnalytic; + +class RooMinimizerFcnSemiAnalytic : + //public ROOT::Math::MultiNumGradFunction { --> need to be constructed on a function. With the nominal Fcn? + public ROOT::Math::IMultiGradFunction { + //public ROOT::Math::IBaseFunctionMultiDim { + + public: + + RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic *context, const std::map& knownDerivatives, bool verbose = false); + RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other); + virtual ~RooMinimizerFcnSemiAnalytic(); + + ROOT::Math::IBaseFunctionMultiDim* Clone() const override; + unsigned int NDim() const override{ return _nDim; } + + RooArgList* GetFloatParamList() { return _floatParamList; } + RooArgList* GetConstParamList() { return _constParamList; } + RooArgList* GetInitFloatParamList() { return _initFloatParamList; } + RooArgList* GetInitConstParamList() { return _initConstParamList; } + + void SetEvalErrorWall(Bool_t flag) { _doEvalErrorWall = flag ; } + void SetPrintEvalErrors(Int_t numEvalErrors) { _printEvalErrors = numEvalErrors ; } + Bool_t SetLogFile(const char* inLogfile); + std::ofstream* GetLogFile() { return _logfile; } + void SetVerbose(Bool_t flag=kTRUE) { _verbose = flag ; } + + Double_t& GetMaxFCN() { return _maxFCN; } + Int_t GetNumInvalidNLL() { return _numBadNLL; } + + Bool_t Synchronize(std::vector& parameters, + Bool_t optConst, Bool_t verbose); + void BackProp(const ROOT::Fit::FitResult &results); + void ApplyCovarianceMatrix(TMatrixDSym& V); + + Int_t evalCounter() const { return _evalCounter ; } + void zeroEvalCount() { _evalCounter = 0 ; } + + + private: + + Double_t GetPdfParamVal(Int_t index); + Double_t GetPdfParamErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t value); + void ClearPdfParamAsymErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal); + + inline Bool_t SetPdfParamVal(const Int_t &index, const Double_t &value) const; + + + double DoEval(const double * x) const override; + double DoDerivative(const double * x,unsigned int icoord) const override; + // semi-analytic. The following function is called by DoDerivative if analytical derivative is not present. + virtual double DoNumericalDerivative(const double * x,int icoord) const; + void updateFloatVec() ; + +private: + + mutable Int_t _evalCounter ; + + RooAbsReal *_funct; + // for name -> Derivative + std::map _knownDerivatives; // does not own them. Otherwise needs to design ad hoc copy constructors. + + RooMinimizerSemiAnalytic *_context; + + mutable double _maxFCN; + mutable int _numBadNLL; + mutable int _printEvalErrors; + Bool_t _doEvalErrorWall; + + int _nDim; + std::ofstream *_logfile; + bool _verbose; + + RooArgList* _floatParamList; + std::vector _floatParamVec ; + std::vector _derivParamVec ; // vector of derivative functions - algined with the one above + // + RooArgList* _constParamList; + RooArgList* _initFloatParamList; + RooArgList* _initConstParamList; + + int _useNumDerivatives{1};// 0 = ROOT (not-impl), 1 re-implementation of gsl + +}; + +#endif +#endif diff --git a/interface/RooMinimizerSemiAnalytic.h b/interface/RooMinimizerSemiAnalytic.h new file mode 100644 index 00000000000..2ed1857d3c5 --- /dev/null +++ b/interface/RooMinimizerSemiAnalytic.h @@ -0,0 +1,140 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + + +/* ~Copy Of RooMINIMIZER, + * the difference is that is templated in RooMinimizerFcn. i + * An instance of RooMinimizerCopy is present as typedef for compatibility. + */ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifndef ROO_MINIMIZER_SEMIANALYTIC +#define ROO_MINIMIZER_SEMIANALYTIC + +#include "TObject.h" +#include "TStopwatch.h" +#include +#include "TMatrixDSymfwd.h" + + +#include "Fit/Fitter.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerFcnSemiAnalytic.h" +#include "RooMinimizerFcn.h" + +class RooAbsReal ; +class RooFitResult ; +class RooArgList ; +class RooRealVar ; +class RooArgSet ; +class TH2F ; +class RooPlot ; + +class RooMinimizerSemiAnalytic : public TObject { +public: + + RooMinimizerSemiAnalytic(RooAbsReal& function,const std::map& knownDerivatives) ; + virtual ~RooMinimizerSemiAnalytic() ; + + enum Strategy { Speed=0, Balance=1, Robustness=2 } ; + enum PrintLevel { None=-1, Reduced=0, Normal=1, ExtraForProblem=2, Maximum=3 } ; + void setStrategy(Int_t strat) ; + void setErrorLevel(Double_t level) ; + void setEps(Double_t eps) ; + void optimizeConst(Int_t flag) ; + void setEvalErrorWall(Bool_t flag) { fitterFcn()->SetEvalErrorWall(flag); } + void setOffsetting(Bool_t flag) ; + void setMaxIterations(Int_t n) ; + void setMaxFunctionCalls(Int_t n) ; + + RooFitResult* fit(const char* options) ; + + Int_t migrad() ; + Int_t hesse() ; + Int_t minos() ; + Int_t minos(const RooArgSet& minosParamList) ; + Int_t seek() ; + Int_t simplex() ; + Int_t improve() ; + + Int_t minimize(const char* type, const char* alg=0) ; + + RooFitResult* save(const char* name=0, const char* title=0) ; + RooPlot* contour(RooRealVar& var1, RooRealVar& var2, + Double_t n1=1, Double_t n2=2, Double_t n3=0, + Double_t n4=0, Double_t n5=0, Double_t n6=0, unsigned int npoints = 50) ; + + Int_t setPrintLevel(Int_t newLevel) ; + void setPrintEvalErrors(Int_t numEvalErrors) { fitterFcn()->SetPrintEvalErrors(numEvalErrors); } + void setVerbose(Bool_t flag=kTRUE) { _verbose = flag ; fitterFcn()->SetVerbose(flag); } + void setProfile(Bool_t flag=kTRUE) { _profile = flag ; } + Bool_t setLogFile(const char* logf=0) { return fitterFcn()->SetLogFile(logf); } + + void setMinimizerType(const char* type) ; + + static void cleanup() ; + static RooFitResult* lastMinuitFit(const RooArgList& varList=RooArgList()) ; + + void saveStatus(const char* label, Int_t status) { _statusHistory.push_back(std::pair(label,status)) ; } + + Int_t evalCounter() const { return fitterFcn()->evalCounter() ; } + void zeroEvalCount() { fitterFcn()->zeroEvalCount() ; } + + ROOT::Fit::Fitter* fitter() ; + const ROOT::Fit::Fitter* fitter() const ; + +protected: + + friend class RooAbsPdf ; + void applyCovarianceMatrix(TMatrixDSym& V) ; + + void profileStart() ; + void profileStop() ; + + inline Int_t getNPar() const { return fitterFcn()->NDim() ; } + inline std::ofstream* logfile() { return fitterFcn()->GetLogFile(); } + inline Double_t& maxFCN() { return fitterFcn()->GetMaxFCN() ; } + + const RooMinimizerFcnSemiAnalytic* fitterFcn() const { return ( fitter()->GetFCN() ? (dynamic_cast ( fitter()->GetFCN()) ): dynamic_cast(_fcn) ) ; } + RooMinimizerFcnSemiAnalytic* fitterFcn() { return ( fitter()->GetFCN() ? (dynamic_cast (fitter()->GetFCN())) : dynamic_cast(_fcn) ) ; } + +private: + + Int_t _printLevel ; + Int_t _status ; + Bool_t _optConst ; + Bool_t _profile ; + RooAbsReal* _func ; + + Bool_t _verbose ; + TStopwatch _timer ; + TStopwatch _cumulTimer ; + Bool_t _profileStart ; + + TMatrixDSym* _extV ; + + RooMinimizerFcnSemiAnalytic *_fcn; // Here it is the difference + std::string _minimizerType; + + static ROOT::Fit::Fitter *_theFitter ; + + std::vector > _statusHistory ; + + RooMinimizerSemiAnalytic(const RooMinimizerSemiAnalytic&) ; + +public: + ClassDef(RooMinimizerSemiAnalytic,0) // RooMinimizerSemiAnalytic this should be a title? +} ; + + +#endif + +#endif diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx new file mode 100644 index 00000000000..62e9082598c --- /dev/null +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -0,0 +1,716 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on roofit code. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +////////////////////////////////////////////////////////////////////////////// +// +// RooMinimizerFcnSemiAnalytic is an interface class to the ROOT::Math function +// for minization. +// + +#include + +#include "RooFit.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerFcnSemiAnalytic.h" + +#include "Riostream.h" + +#include "TIterator.h" +#include "TClass.h" + +#include "RooAbsArg.h" +#include "RooAbsPdf.h" +#include "RooArgSet.h" +#include "RooRealVar.h" +#include "RooAbsRealLValue.h" +#include "RooMsgService.h" + +#include "RooMinimizer.h" +#include + +using namespace std; + +RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic* context, const std::map& knownDerivatives, bool verbose) : + _funct(funct), + _knownDerivatives (knownDerivatives), + _context(context), + // Reset the *largest* negative log-likelihood value we have seen so far + _maxFCN(-1e30), _numBadNLL(0), + _printEvalErrors(10), _doEvalErrorWall(kTRUE), + _nDim(0), _logfile(0), + _verbose(verbose) +{ + + _evalCounter = 0 ; + + // Examine parameter list + RooArgSet* paramSet = _funct->getParameters(RooArgSet()); + RooArgList paramList(*paramSet); + delete paramSet; + + _floatParamList = (RooArgList*) paramList.selectByAttrib("Constant",kFALSE); + if (_floatParamList->getSize()>1) { + _floatParamList->sort(); + } + _floatParamList->setName("floatParamList"); + + _constParamList = (RooArgList*) paramList.selectByAttrib("Constant",kTRUE); + if (_constParamList->getSize()>1) { + _constParamList->sort(); + } + _constParamList->setName("constParamList"); + + // Remove all non-RooRealVar parameters from list (MINUIT cannot handle them) + TIterator* pIter = _floatParamList->createIterator(); + RooAbsArg* arg; + while ((arg=(RooAbsArg*)pIter->Next())) { + if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic: removing parameter " + << arg->GetName() + << " from list because it is not of type RooRealVar" << endl; + _floatParamList->remove(*arg); + } + } + delete pIter; + + _nDim = _floatParamList->getSize(); + + updateFloatVec() ; + + // Save snapshot of initial lists + _initFloatParamList = (RooArgList*) _floatParamList->snapshot(kFALSE) ; + _initConstParamList = (RooArgList*) _constParamList->snapshot(kFALSE) ; + +} + + + +RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other) : ROOT::Math::IMultiGradFunction(other), + _evalCounter(other._evalCounter), + _funct(other._funct), + _knownDerivatives(other._knownDerivatives), + _context(other._context), + _maxFCN(other._maxFCN), + _numBadNLL(other._numBadNLL), + _printEvalErrors(other._printEvalErrors), + _doEvalErrorWall(other._doEvalErrorWall), + _nDim(other._nDim), + _logfile(other._logfile), + _verbose(other._verbose), + _floatParamVec(other._floatParamVec), + _derivParamVec(other._derivParamVec) +{ + _floatParamList = new RooArgList(*other._floatParamList) ; + _constParamList = new RooArgList(*other._constParamList) ; + _initFloatParamList = (RooArgList*) other._initFloatParamList->snapshot(kFALSE) ; + _initConstParamList = (RooArgList*) other._initConstParamList->snapshot(kFALSE) ; +} + + +RooMinimizerFcnSemiAnalytic::~RooMinimizerFcnSemiAnalytic() +{ + delete _floatParamList; + delete _initFloatParamList; + delete _constParamList; + delete _initConstParamList; +} + + +ROOT::Math::IBaseFunctionMultiDim* RooMinimizerFcnSemiAnalytic::Clone() const +{ + return new RooMinimizerFcnSemiAnalytic(*this) ; +} + + +Bool_t RooMinimizerFcnSemiAnalytic::Synchronize(std::vector& parameters, + Bool_t optConst, Bool_t verbose) +{ + + // Internal function to synchronize TMinimizer with current + // information in RooAbsReal function parameters + + Bool_t constValChange(kFALSE) ; + Bool_t constStatChange(kFALSE) ; + + Int_t index(0) ; + + // Handle eventual migrations from constParamList -> floatParamList + for(index= 0; index < _constParamList->getSize() ; index++) { + + RooRealVar *par= dynamic_cast(_constParamList->at(index)) ; + if (!par) continue ; + + RooRealVar *oldpar= dynamic_cast(_initConstParamList->at(index)) ; + if (!oldpar) continue ; + + // Test if constness changed + if (!par->isConstant()) { + + // Remove from constList, add to floatList + _constParamList->remove(*par) ; + _floatParamList->add(*par) ; + _initFloatParamList->addClone(*oldpar) ; + _initConstParamList->remove(*oldpar) ; + constStatChange=kTRUE ; + _nDim++ ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now floating." << endl ; + } + } + + // Test if value changed + if (par->getVal()!= oldpar->getVal()) { + constValChange=kTRUE ; + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of constant parameter " + << par->GetName() + << " changed from " << oldpar->getVal() << " to " + << par->getVal() << endl ; + } + } + + } + + // Update reference list + *_initConstParamList = *_constParamList ; + + // Synchronize MINUIT with function state + // Handle floatParamList + for(index= 0; index < _floatParamList->getSize(); index++) { + RooRealVar *par= dynamic_cast(_floatParamList->at(index)) ; + + if (!par) continue ; + + Double_t pstep(0) ; + Double_t pmin(0) ; + Double_t pmax(0) ; + + if(!par->isConstant()) { + + // Verify that floating parameter is indeed of type RooRealVar + if (!par->IsA()->InheritsFrom(RooRealVar::Class())) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::fit: Error, non-constant parameter " + << par->GetName() + << " is not of type RooRealVar, skipping" << endl ; + _floatParamList->remove(*par); + index--; + _nDim--; + continue ; + } + + // Set the limits, if not infinite + if (par->hasMin() ) + pmin = par->getMin(); + if (par->hasMax() ) + pmax = par->getMax(); + + // Calculate step size + pstep = par->getError(); + if(pstep <= 0) { + // Floating parameter without error estitimate + if (par->hasMin() && par->hasMax()) { + pstep= 0.1*(pmax-pmin); + + // Trim default choice of error if within 2 sigma of limit + if (pmax - par->getVal() < 2*pstep) { + pstep = (pmax - par->getVal())/2 ; + } else if (par->getVal() - pmin < 2*pstep) { + pstep = (par->getVal() - pmin )/2 ; + } + + // If trimming results in zero error, restore default + if (pstep==0) { + pstep= 0.1*(pmax-pmin); + } + + } else { + pstep=1 ; + } + if(verbose) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::synchronize: WARNING: no initial error estimate available for " + << par->GetName() << ": using " << pstep << endl; + } + } + } else { + pmin = par->getVal() ; + pmax = par->getVal() ; + } + + // new parameter + if (index>=Int_t(parameters.size())) { + + if (par->hasMin() && par->hasMax()) { + parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(), + par->getVal(), + pstep, + pmin,pmax)); + } + else { + parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(), + par->getVal(), + pstep)); + if (par->hasMin() ) + parameters.back().SetLowerLimit(pmin); + else if (par->hasMax() ) + parameters.back().SetUpperLimit(pmax); + } + + continue; + + } + + Bool_t oldFixed = parameters[index].IsFixed(); + Double_t oldVar = parameters[index].Value(); + Double_t oldVerr = parameters[index].StepSize(); + Double_t oldVlo = parameters[index].LowerLimit(); + Double_t oldVhi = parameters[index].UpperLimit(); + + if (par->isConstant() && !oldFixed) { + + // Parameter changes floating -> constant : update only value if necessary + if (oldVar!=par->getVal()) { + parameters[index].SetValue(par->getVal()); + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of parameter " + << par->GetName() << " changed from " << oldVar + << " to " << par->getVal() << endl ; + } + } + parameters[index].Fix(); + constStatChange=kTRUE ; + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now fixed." << endl ; + } + + } else if (par->isConstant() && oldFixed) { + + // Parameter changes constant -> constant : update only value if necessary + if (oldVar!=par->getVal()) { + parameters[index].SetValue(par->getVal()); + constValChange=kTRUE ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of fixed parameter " + << par->GetName() << " changed from " << oldVar + << " to " << par->getVal() << endl ; + } + } + + } else { + // Parameter changes constant -> floating + if (!par->isConstant() && oldFixed) { + parameters[index].Release(); + constStatChange=kTRUE ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now floating." << endl ; + } + } + + // Parameter changes constant -> floating : update all if necessary + if (oldVar!=par->getVal() || oldVlo!=pmin || oldVhi != pmax || oldVerr!=pstep) { + parameters[index].SetValue(par->getVal()); + parameters[index].SetStepSize(pstep); + if (par->hasMin() && par->hasMax() ) + parameters[index].SetLimits(pmin,pmax); + else if (par->hasMin() ) + parameters[index].SetLowerLimit(pmin); + else if (par->hasMax() ) + parameters[index].SetUpperLimit(pmax); + } + + // Inform user about changes in verbose mode + if (verbose) { + // if ierr<0, par was moved from the const list and a message was already printed + + if (oldVar!=par->getVal()) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of parameter " + << par->GetName() << " changed from " << oldVar << " to " + << par->getVal() << endl ; + } + if (oldVlo!=pmin || oldVhi!=pmax) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: limits of parameter " + << par->GetName() << " changed from [" << oldVlo << "," << oldVhi + << "] to [" << pmin << "," << pmax << "]" << endl ; + } + + // If oldVerr=0, then parameter was previously fixed + if (oldVerr!=pstep && oldVerr!=0) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: error/step size of parameter " + << par->GetName() << " changed from " << oldVerr << " to " << pstep << endl ; + } + } + } + } + + if (optConst) { + if (constStatChange) { + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: set of constant parameters changed, rerunning const optimizer" << endl ; + _funct->constOptimizeTestStatistic(RooAbsArg::ConfigChange) ; + } else if (constValChange) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: constant parameter values changed, rerunning const optimizer" << endl ; + _funct->constOptimizeTestStatistic(RooAbsArg::ValueChange) ; + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + + } + + updateFloatVec() ; + + return 0 ; + +} + +Double_t RooMinimizerFcnSemiAnalytic::GetPdfParamVal(Int_t index) +{ + // Access PDF parameter value by ordinal index (needed by MINUIT) + + return ((RooRealVar*)_floatParamList->at(index))->getVal() ; +} + +Double_t RooMinimizerFcnSemiAnalytic::GetPdfParamErr(Int_t index) +{ + // Access PDF parameter error by ordinal index (needed by MINUIT) + return ((RooRealVar*)_floatParamList->at(index))->getError() ; +} + + +void RooMinimizerFcnSemiAnalytic::SetPdfParamErr(Int_t index, Double_t value) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->setError(value) ; +} + + + +void RooMinimizerFcnSemiAnalytic::ClearPdfParamAsymErr(Int_t index) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->removeAsymError() ; +} + + +void RooMinimizerFcnSemiAnalytic::SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->setAsymError(loVal,hiVal) ; +} + + +void RooMinimizerFcnSemiAnalytic::BackProp(const ROOT::Fit::FitResult &results) +{ + // Transfer MINUIT fit results back into RooFit objects + + for (Int_t index= 0; index < _nDim; index++) { + Double_t value = results.Value(index); + SetPdfParamVal(index, value); + + // Set the parabolic error + Double_t err = results.Error(index); + SetPdfParamErr(index, err); + + Double_t eminus = results.LowerError(index); + Double_t eplus = results.UpperError(index); + + if(eplus > 0 || eminus < 0) { + // Store the asymmetric error, if it is available + SetPdfParamErr(index, eminus,eplus); + } else { + // Clear the asymmetric error + ClearPdfParamAsymErr(index) ; + } + } + +} + +Bool_t RooMinimizerFcnSemiAnalytic::SetLogFile(const char* inLogfile) +{ + // Change the file name for logging of a RooMinimizer of all MINUIT steppings + // through the parameter space. If inLogfile is null, the current log file + // is closed and logging is stopped. + + if (_logfile) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::setLogFile: closing previous log file" << endl ; + _logfile->close() ; + delete _logfile ; + _logfile = 0 ; + } + _logfile = new ofstream(inLogfile) ; + if (!_logfile->good()) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::setLogFile: cannot open file " << inLogfile << endl ; + _logfile->close() ; + delete _logfile ; + _logfile= 0; + } + + return kFALSE ; + +} + + +void RooMinimizerFcnSemiAnalytic::ApplyCovarianceMatrix(TMatrixDSym& V) +{ + // Apply results of given external covariance matrix. i.e. propagate its errors + // to all RRV parameter representations and give this matrix instead of the + // HESSE matrix at the next save() call + + for (Int_t i=0 ; i<_nDim ; i++) { + // Skip fixed parameters + if (_floatParamList->at(i)->isConstant()) { + continue ; + } + SetPdfParamErr(i, sqrt(V(i,i))) ; + } + +} + + +Bool_t RooMinimizerFcnSemiAnalytic::SetPdfParamVal(const Int_t &index, const Double_t &value) const +{ + //RooRealVar* par = (RooRealVar*)_floatParamList->at(index); + RooRealVar* par = (RooRealVar*)_floatParamVec[index] ; + + if (par->getVal()!=value) { + if (_verbose) cout << par->GetName() << "=" << value << ", " ; + + par->setVal(value); + return kTRUE; + } + + return kFALSE; +} + + + +//////////////////////////////////////////////////////////////////////////////// + +void RooMinimizerFcnSemiAnalytic::updateFloatVec() +{ + // derivative vector and float vectors needs to be aligned + _floatParamVec.clear() ; + _derivParamVec.clear(); + + RooFIter iter = _floatParamList->fwdIterator() ; + RooAbsArg* arg ; + _floatParamVec = std::vector(_floatParamList->getSize()) ; + _derivParamVec = std::vector(_floatParamList->getSize()) ; + Int_t i(0) ; + while((arg=iter.next())) { + auto derivative=_knownDerivatives.find(arg->GetName()); + if (derivative == _knownDerivatives.end()){ + if(_verbose) std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "derivative for param "<GetName()<<" not found. I will use numerical derivatives."<GetName()<<" found with name "<second->GetName()<<"."<second; + } + _floatParamVec[i++] = arg ; + } +} + + + +double RooMinimizerFcnSemiAnalytic::DoEval(const double *x) const +{ + + // Set the parameter values for this iteration + for (int index = 0; index < _nDim; index++) { + if (_logfile) (*_logfile) << x[index] << " " ; + SetPdfParamVal(index,x[index]); + } + + // Calculate the function for these parameters + RooAbsReal::setHideOffset(kFALSE) ; + double fvalue = _funct->getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + + if (RooAbsReal::numEvalErrors()>0 || fvalue > 1e30) { + + if (_printEvalErrors>=0) { + + if (_doEvalErrorWall) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic: Minimized function has error status." << endl + << "Returning maximum FCN so far (" << _maxFCN + << ") to force MIGRAD to back out of this region. Error log follows" << endl ; + } else { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic: Minimized function has error status but is ignored" << endl ; + } + + Bool_t first(kTRUE) ; + /*ooccoutW(_context,Minimization)*/std::cout << "Parameter values: " ; + //for (const auto par : *_floatParamList) { //new loop + for (int ip =0 ;ip < _floatParamList->getSize();++ip){ + const auto par=_floatParamList->at(ip); + auto var = static_cast(par); + if (first) { first = kFALSE ; } else /*ooccoutW(_context,Minimization)*/std::cout << ", " ; + /*ooccoutW(_context,Minimization)*/std::cout << var->GetName() << "=" << var->getVal() ; + } + /*ooccoutW(_context,Minimization)*/std::cout << endl ; + + RooAbsReal::printEvalErrors(/*ooccoutW(_context,Minimization)*/std::cout,_printEvalErrors) ; + /*ooccoutW(_context,Minimization)*/std::cout << endl ; + } + + if (_doEvalErrorWall) { + fvalue = _maxFCN+1; + } + + RooAbsReal::clearEvalErrorLog() ; + _numBadNLL++ ; + } else { + _maxFCN = std::max(fvalue, _maxFCN); + } + + // Optional logging + if (_logfile) + (*_logfile) << setprecision(15) << fvalue << setprecision(4) << endl; + if (_verbose) { + cout << "\nprevFCN" << (_funct->isOffsetting()?"-offset":"") << " = " << setprecision(10) + << fvalue << setprecision(4) << " " ; + cout.flush() ; + } + + _evalCounter++ ; + + return fvalue; +} + +double RooMinimizerFcnSemiAnalytic::DoDerivative(const double * x, unsigned int icoord) const{ + + // Set the parameter values for this iteration + for (int index = 0; index < _nDim; index++) { + if (_logfile) (*_logfile) << x[index] << " " ; + SetPdfParamVal(index,x[index]); + } + + // compute Derivative + if (_derivParamVec.size()getVal(); + +} + +double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int icoord) const +{ + // ROOT + if (_useNumDerivatives==0) { + //return ROOT::Math::MultiNumGradFunction::DoDerivative(x,icoord); + //copied from ROOT::Math::MultiNumGradFunction::DoDerivative + + // calculate derivative using mathcore derivator class + // step size can be changes using SetDerivPrecision() + // this also use gsl. + + //static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + //double x0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + //double step = (x0>0) ? std::max( fgEps* x0, 8.0*kPrecision*(x0 + kPrecision) ) : kPrecision; + //return ROOT::Math::Derivator::Eval(*fFunc, x, icoord, step); + throw std::logic_error("Unimplemented"); + } + if (_useNumDerivatives==1){ + /* + * + ROOT use this has step above, maybe copy it dynamically? TODO + double MultiNumGradFunction::fgEps = 0.001; + + double MultiNumGradFunction::DoDerivative (const double * x, unsigned int icoord ) const { + // calculate derivative using mathcore derivator class + // step size can be changes using SetDerivPrecision() + + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + double x0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double step = (x0>0) ? std::max( fgEps* x0, 8.0*kPrecision*(x0 + kPrecision) ) : kPrecision; + */ + + // from gsl_deriv_central -> central_deriv + /* Compute the derivative using the 5-point rule (x-h, x-h/2, x, + x+h/2, x+h). Note that the central point is not used. + + Compute the error using the difference between the 5-point and + the 3-point rule (x-h,x,x+h). Again the central point is not + used. */ + + // smallest number representable in doubles + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + static double fgEps = 0.001; + + double ax0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double h = (ax0>0) ? std::max( fgEps* ax0, 8.0*kPrecision*(ax0 + kPrecision) ) : kPrecision; + //static const double h=std::numeric_limits::epsilon()*8.; + const double hhalf=h/2.; + + //double f0 = _funct->getVal(); + //double x0 = x[icoord]; + + SetPdfParamVal(icoord,x[icoord]-h); + double fm1 = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+h); + double fp1 = _funct->getVal(); + + SetPdfParamVal(icoord,x[icoord]-hhalf); + double fmh = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+hhalf); + double fph = _funct->getVal(); + + double r3 = 0.5 * (fp1 - fm1); + double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3; + + // don't do counts I'm not using. + //double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON; + //double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3; + + /* The next term is due to finite precision in x+h = O (eps * x) */ + + //double dy = std::max (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON; + + /* The truncation error in the r5 approximation itself is O(h^4). + However, for safety, we estimate the error from r5-r3, which is + O(h^2). By scaling h we will minimise this estimated error, not + the actual truncation error in r5. */ + + //*result = r5 / h; + //*abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */ + //*abserr_round = fabs (e5 / h) + dy; /* Rounding error (cancellations) */ + return r5 /h ; + } + throw std::runtime_error(Form("Numerical derivatives method un-implemented: %d",_useNumDerivatives)); +} + + +#endif + diff --git a/src/RooMinimizerSemiAnalytic.cc b/src/RooMinimizerSemiAnalytic.cc new file mode 100644 index 00000000000..08f2dae6224 --- /dev/null +++ b/src/RooMinimizerSemiAnalytic.cc @@ -0,0 +1,942 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +/** +\file RooMinimizerSemiAnalytic.cxx +\class RooMinimizerSemiAnalytic +\ingroup Roofitcore + +RooMinimizerSemiAnalytic is a copy of RooMinimizer. +RooMinimizerSemiAnalytic is templated in RooMinimizerFcn. +An instance is present as typedef definined on 'RooMinimizerFcn', that should act as an exact copy of RooMinimizer. +**/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifdef ROO_FIT_RESULT + #error I need to import RooFitResult +#else + #include "HiggsAnalysis/CombinedLimit/interface/RooFitResult.h" +#endif + +#include "RooFit.h" +#include "Riostream.h" + +#include "TClass.h" + +#include + +#include "TH1.h" +#include "TH2.h" +#include "TMarker.h" +#include "TGraph.h" +#include "Fit/FitConfig.h" +#include "TStopwatch.h" +#include "TDirectory.h" +#include "TMatrixDSym.h" + +#include "RooArgSet.h" +#include "RooArgList.h" +#include "RooAbsReal.h" +#include "RooAbsRealLValue.h" +#include "RooRealVar.h" +#include "RooAbsPdf.h" +#include "RooSentinel.h" +#include "RooMsgService.h" +#include "RooPlot.h" + + +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include "RooFitResult.h" + +#include "Math/Minimizer.h" + +#if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3) +char* operator+( streampos&, char* ); +#endif + +using namespace std; + +//ClassImp(RooMinimizerSemiAnalytic); + + +ROOT::Fit::Fitter *RooMinimizerSemiAnalytic::_theFitter = 0 ; +//ROOT::Fit::Fitter *RooMinimizerSemiAnalytic::_theFitter = 0 ; + + + +//////////////////////////////////////////////////////////////////////////////// +/// Cleanup method called by atexit handler installed by RooSentinel +/// to delete all global heap objects when the program is terminated + +void RooMinimizerSemiAnalytic::cleanup() +{ + if (_theFitter) { + delete _theFitter ; + _theFitter =0 ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Construct MINUIT interface to given function. Function can be anything, +/// but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2 +/// (implemented by RooChi2Var). Other frequent use cases are a RooAddition +/// of a RooNLLVar plus a penalty or constraint term. This class propagates +/// all RooFit information (floating parameters, their values and errors) +/// to MINUIT before each MINUIT call and propagates all MINUIT information +/// back to the RooFit object at the end of each call (updated parameter +/// values, their (asymmetric errors) etc. The default MINUIT error level +/// for HESSE and MINOS error analysis is taken from the defaultErrorLevel() +/// value of the input function. + +RooMinimizerSemiAnalytic::RooMinimizerSemiAnalytic(RooAbsReal& function,const std::map& knownDerivatives) +{ + RooSentinel::activate() ; + + // Store function reference + _extV = 0 ; + _func = &function ; + _optConst = kFALSE ; + _verbose = kFALSE ; + _profile = kFALSE ; + _profileStart = kFALSE ; + _printLevel = 1 ; + _minimizerType = "Minuit"; // default minimizer + + if (_theFitter) delete _theFitter ; + _theFitter = new ROOT::Fit::Fitter; + _fcn = new RooMinimizerFcnSemiAnalytic(_func,this,knownDerivatives,_verbose); + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + setEps(1.0); // default tolerance + // default max number of calls + _theFitter->Config().MinimizerOptions().SetMaxIterations(500*_fcn->NDim()); + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500*_fcn->NDim()); + + // Shut up for now + setPrintLevel(-1) ; + + // Use +0.5 for 1-sigma errors + setErrorLevel(_func->defaultErrorLevel()) ; + + // Declare our parameters to MINUIT + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + + // Now set default verbosity + if (RooMsgService::instance().silentMode()) { + setPrintLevel(-1) ; + } else { + setPrintLevel(1) ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Destructor + +RooMinimizerSemiAnalytic::~RooMinimizerSemiAnalytic() +{ + if (_extV) { + delete _extV ; + } + + if (_fcn) { + delete _fcn; + } + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change MINUIT strategy to istrat. Accepted codes +/// are 0,1,2 and represent MINUIT strategies for dealing +/// most efficiently with fast FCNs (0), expensive FCNs (2) +/// and 'intermediate' FCNs (1) + +void RooMinimizerSemiAnalytic::setStrategy(Int_t istrat) +{ + _theFitter->Config().MinimizerOptions().SetStrategy(istrat); + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change maximum number of MINUIT iterations +/// (RooMinimizerSemiAnalytic default 500 * #parameters) + +void RooMinimizerSemiAnalytic::setMaxIterations(Int_t n) +{ + _theFitter->Config().MinimizerOptions().SetMaxIterations(n); +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change maximum number of likelihood function calss from MINUIT +/// (RooMinimizerSemiAnalytic default 500 * #parameters) + +void RooMinimizerSemiAnalytic::setMaxFunctionCalls(Int_t n) +{ + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(n); +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Set the level for MINUIT error analysis to the given +/// value. This function overrides the default value +/// that is taken in the RooMinimizerSemiAnalytic constructor from +/// the defaultErrorLevel() method of the input function + +void RooMinimizerSemiAnalytic::setErrorLevel(Double_t level) +{ + _theFitter->Config().MinimizerOptions().SetErrorDef(level); + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change MINUIT epsilon + +void RooMinimizerSemiAnalytic::setEps(Double_t eps) +{ + _theFitter->Config().MinimizerOptions().SetTolerance(eps); + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Enable internal likelihood offsetting for enhanced numeric precision + +void RooMinimizerSemiAnalytic::setOffsetting(Bool_t flag) +{ + _func->enableOffsetting(flag) ; +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Choose the minimzer algorithm. + +void RooMinimizerSemiAnalytic::setMinimizerType(const char* type) +{ + _minimizerType = type; +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Return underlying ROOT fitter object + +ROOT::Fit::Fitter* RooMinimizerSemiAnalytic::fitter() +{ + return _theFitter ; +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Return underlying ROOT fitter object + +const ROOT::Fit::Fitter* RooMinimizerSemiAnalytic::fitter() const +{ + return _theFitter ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Parse traditional RooAbsPdf::fitTo driver options +/// +/// m - Run Migrad only +/// h - Run Hesse to estimate errors +/// v - Verbose mode +/// l - Log parameters after each Minuit steps to file +/// t - Activate profile timer +/// r - Save fit result +/// 0 - Run Migrad with strategy 0 + +RooFitResult* RooMinimizerSemiAnalytic::fit(const char* options) +{ + TString opts(options) ; + opts.ToLower() ; + + // Initial configuration + if (opts.Contains("v")) setVerbose(1) ; + if (opts.Contains("t")) setProfile(1) ; + if (opts.Contains("l")) setLogFile(Form("%s.log",_func->GetName())) ; + if (opts.Contains("c")) optimizeConst(1) ; + + // Fitting steps + if (opts.Contains("0")) setStrategy(0) ; + migrad() ; + if (opts.Contains("0")) setStrategy(1) ; + if (opts.Contains("h")||!opts.Contains("m")) hesse() ; + if (!opts.Contains("m")) minos() ; + + return (opts.Contains("r")) ? save() : 0 ; +} + + + + +//////////////////////////////////////////////////////////////////////////////// + +Int_t RooMinimizerSemiAnalytic::minimize(const char* type, const char* alg) +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + + _theFitter->Config().SetMinimizer(type,alg); + + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + bool ret = _theFitter->FitFCN(*_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINIMIZE",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MIGRAD. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::migrad() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad"); + bool ret = _theFitter->FitFCN(*_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MIGRAD",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute HESSE. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::hesse() +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::hesse: Error, run Migrad before Hesse!" + << endl ; + _status = -1; + } + else { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateHessErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("HESSE",_status) ; + + } + + return _status ; + +} + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MINOS. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::minos() +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::minos: Error, run Migrad before Minos!" + << endl ; + _status = -1; + } + else { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateMinosErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINOS",_status) ; + + } + + return _status ; + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MINOS for given list of parameters. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::minos(const RooArgSet& minosParamList) +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::minos: Error, run Migrad before Minos!" + << endl ; + _status = -1; + } + else if (minosParamList.getSize()>0) { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + // get list of parameters for Minos + TIterator* aIter = minosParamList.createIterator() ; + RooAbsArg* arg ; + std::vector paramInd; + while((arg=(RooAbsArg*)aIter->Next())) { + RooAbsArg* par = _fcn->GetFloatParamList()->find(arg->GetName()); + if (par && !par->isConstant()) { + Int_t index = _fcn->GetFloatParamList()->index(par); + paramInd.push_back(index); + } + } + delete aIter ; + + if (paramInd.size()) { + // set the parameter indeces + _theFitter->Config().SetMinosErrors(paramInd); + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateMinosErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + // to avoid that following minimization computes automatically the Minos errors + _theFitter->Config().SetMinosErrors(kFALSE); + + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINOS",_status) ; + + } + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute SEEK. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::seek() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek"); + bool ret = _theFitter->FitFCN(*_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("SEEK",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute SIMPLEX. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::simplex() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex"); + bool ret = _theFitter->FitFCN(*_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("SEEK",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute IMPROVE. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::improve() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved"); + bool ret = _theFitter->FitFCN(*_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("IMPROVE",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change the MINUIT internal printing level + +Int_t RooMinimizerSemiAnalytic::setPrintLevel(Int_t newLevel) +{ + Int_t ret = _printLevel ; + _theFitter->Config().MinimizerOptions().SetPrintLevel(newLevel+1); + _printLevel = newLevel+1 ; + return ret ; +} + +//////////////////////////////////////////////////////////////////////////////// +/// If flag is true, perform constant term optimization on +/// function being minimized. + +void RooMinimizerSemiAnalytic::optimizeConst(Int_t flag) +{ + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + + if (_optConst && !flag){ + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: deactivating const optimization" << endl ; + _func->constOptimizeTestStatistic(RooAbsArg::DeActivate) ; + _optConst = flag ; + } else if (!_optConst && flag) { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: activating const optimization" << endl ; + _func->constOptimizeTestStatistic(RooAbsArg::Activate,flag>1) ; + _optConst = flag ; + } else if (_optConst && flag) { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: const optimization already active" << endl ; + } else { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: const optimization wasn't active" << endl ; + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Save and return a RooFitResult snaphot of current minimizer status. +/// This snapshot contains the values of all constant parameters, +/// the value of all floating parameters at RooMinimizerSemiAnalytic construction and +/// after the last MINUIT operation, the MINUIT status, variance quality, +/// EDM setting, number of calls with evaluation problems, the minimized +/// function value and the full correlation matrix + +RooFitResult* RooMinimizerSemiAnalytic::save(const char* userName, const char* userTitle) +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::save: Error, run minimization before!" + << endl ; + return 0; + } + + TString name,title ; + name = userName ? userName : Form("%s", _func->GetName()) ; + title = userTitle ? userTitle : Form("%s", _func->GetTitle()) ; + RooFitResult* fitRes = new RooFitResult(name,title) ; + + // Move eventual fixed parameters in floatList to constList + Int_t i ; + RooArgList saveConstList(*(_fcn->GetConstParamList())) ; + RooArgList saveFloatInitList(*(_fcn->GetInitFloatParamList())) ; + RooArgList saveFloatFinalList(*(_fcn->GetFloatParamList())) ; + for (i=0 ; i<_fcn->GetFloatParamList()->getSize() ; i++) { + RooAbsArg* par = _fcn->GetFloatParamList()->at(i) ; + if (par->isConstant()) { + saveFloatInitList.remove(*saveFloatInitList.find(par->GetName()),kTRUE) ; + saveFloatFinalList.remove(*par) ; + saveConstList.add(*par) ; + } + } + saveConstList.sort() ; + + fitRes->setConstParList(saveConstList) ; + fitRes->setInitParList(saveFloatInitList) ; + + fitRes->setStatus(_status) ; + fitRes->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ; + fitRes->setMinNLL(_theFitter->Result().MinFcnValue()) ; + fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()) ; + fitRes->setEDM(_theFitter->Result().Edm()) ; + fitRes->setFinalParList(saveFloatFinalList) ; + if (!_extV) { + std::vector globalCC; + TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ; + TMatrixDSym covs(_theFitter->Result().Parameters().size()) ; + for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) { + globalCC.push_back(_theFitter->Result().GlobalCC(ic)); + for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) { + corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii); + covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii); + } + } + fitRes->fillCorrMatrix(globalCC,corrs,covs) ; + } else { + fitRes->setCovarianceMatrix(*_extV) ; + } + + fitRes->setStatusHistory(_statusHistory) ; + + return fitRes ; + +} + +//////////////////////////////////////////////////////////////////////////////// +/// Create and draw a TH2 with the error contours in the parameters `var1` and `var2`. +/// \param[in] var1 The first parameter (x axis). +/// \param[in] var2 The second parameter (y axis). +/// \param[in] n1 First contour. +/// \param[in] n2 Optional contour. 0 means don't draw. +/// \param[in] n3 Optional contour. 0 means don't draw. +/// \param[in] n4 Optional contour. 0 means don't draw. +/// \param[in] n5 Optional contour. 0 means don't draw. +/// \param[in] n6 Optional contour. 0 means don't draw. +/// \param[in] npoints Number of points for evaluating the contour. +/// +/// Up to six contours can be drawn using the arguments `n1` to `n6` to request the desired +/// coverage in units of \f$ \sigma = n^2 \cdot \mathrm{ErrorDef} \f$. +/// See ROOT::Math::Minimizer::ErrorDef(). + +RooPlot* RooMinimizerSemiAnalytic::contour(RooRealVar& var1, RooRealVar& var2, + Double_t n1, Double_t n2, Double_t n3, + Double_t n4, Double_t n5, Double_t n6, unsigned int npoints) +{ + + + RooArgList* params = _fcn->GetFloatParamList() ; + RooArgList* paramSave = (RooArgList*) params->snapshot() ; + + // Verify that both variables are floating parameters of PDF + Int_t index1= _fcn->GetFloatParamList()->index(&var1); + if(index1 < 0) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" << GetName() + << ") ERROR: " << var1.GetName() + << " is not a floating parameter of " + << _func->GetName() << endl ; + return 0; + } + + Int_t index2= _fcn->GetFloatParamList()->index(&var2); + if(index2 < 0) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" << GetName() + << ") ERROR: " << var2.GetName() + << " is not a floating parameter of PDF " + << _func->GetName() << endl ; + return 0; + } + + // create and draw a frame + RooPlot* frame = new RooPlot(var1,var2) ; + + // draw a point at the current parameter values + TMarker *point= new TMarker(var1.getVal(), var2.getVal(), 8); + frame->addObject(point) ; + + // check first if a inimizer is available. If not means + // the minimization is not done , so do it + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::contour: Error, run Migrad before contours!" + << endl ; + return frame; + } + + + // remember our original value of ERRDEF + Double_t errdef= _theFitter->GetMinimizer()->ErrorDef(); + + Double_t n[6] ; + n[0] = n1 ; n[1] = n2 ; n[2] = n3 ; n[3] = n4 ; n[4] = n5 ; n[5] = n6 ; + + for (Int_t ic = 0 ; ic<6 ; ic++) { + if(n[ic] > 0) { + + // set the value corresponding to an n1-sigma contour + _theFitter->GetMinimizer()->SetErrorDef(n[ic]*n[ic]*errdef); + + // calculate and draw the contour + Double_t *xcoor = new Double_t[npoints+1]; + Double_t *ycoor = new Double_t[npoints+1]; + bool ret = _theFitter->GetMinimizer()->Contour(index1,index2,npoints,xcoor,ycoor); + + if (!ret) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" + << GetName() + << ") ERROR: MINUIT did not return a contour graph for n=" + << n[ic] << endl ; + } else { + xcoor[npoints] = xcoor[0]; + ycoor[npoints] = ycoor[0]; + TGraph* graph = new TGraph(npoints+1,xcoor,ycoor); + + graph->SetName(Form("contour_%s_n%f",_func->GetName(),n[ic])) ; + graph->SetLineStyle(ic+1) ; + graph->SetLineWidth(2) ; + graph->SetLineColor(kBlue) ; + frame->addObject(graph,"L") ; + } + + delete [] xcoor; + delete [] ycoor; + } + } + + + // restore the original ERRDEF + _theFitter->Config().MinimizerOptions().SetErrorDef(errdef); + + // restore parameter values + *params = *paramSave ; + delete paramSave ; + + return frame ; + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Start profiling timer + +void RooMinimizerSemiAnalytic::profileStart() +{ + if (_profile) { + _timer.Start() ; + _cumulTimer.Start(_profileStart?kFALSE:kTRUE) ; + _profileStart = kTRUE ; + } +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Stop profiling timer and report results of last session + +void RooMinimizerSemiAnalytic::profileStop() +{ + if (_profile) { + _timer.Stop() ; + _cumulTimer.Stop() ; + coutI(Minimization) << "Command timer: " ; _timer.Print() ; + coutI(Minimization) << "Session timer: " ; _cumulTimer.Print() ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Apply results of given external covariance matrix. i.e. propagate its errors +/// to all RRV parameter representations and give this matrix instead of the +/// HESSE matrix at the next save() call + +void RooMinimizerSemiAnalytic::applyCovarianceMatrix(TMatrixDSym& V) +{ + _extV = (TMatrixDSym*) V.Clone() ; + _fcn->ApplyCovarianceMatrix(*_extV); + +} + + +RooFitResult* RooMinimizerSemiAnalytic::lastMinuitFit(const RooArgList& varList) +{ + // Import the results of the last fit performed, interpreting + // the fit parameters as the given varList of parameters. + + if (_theFitter==0 || _theFitter->GetMinimizer()==0) { + oocoutE((TObject*)0,InputArguments) << "RooMinimizerSemiAnalytic::save: Error, run minimization before!" + << endl ; + return 0; + } + + // Verify length of supplied varList + if (varList.getSize()>0 && varList.getSize()!=Int_t(_theFitter->Result().NTotalParameters())) { + oocoutE((TObject*)0,InputArguments) + << "RooMinimizerSemiAnalytic::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl + << " or match the number of variables of the last fit (" + << _theFitter->Result().NTotalParameters() << ")" << endl ; + return 0 ; + } + + + // Verify that all members of varList are of type RooRealVar + TIterator* iter = varList.createIterator() ; + RooAbsArg* arg ; + while((arg=(RooAbsArg*)iter->Next())) { + if (!dynamic_cast(arg)) { + oocoutE((TObject*)0,InputArguments) << "RooMinimizerSemiAnalytic::lastMinuitFit: ERROR: variable '" + << arg->GetName() << "' is not of type RooRealVar" << endl ; + return 0 ; + } + } + delete iter ; + + RooFitResult* res = new RooFitResult("lastMinuitFit","Last MINUIT fit") ; + + // Extract names of fit parameters + // and construct corresponding RooRealVars + RooArgList constPars("constPars") ; + RooArgList floatPars("floatPars") ; + + UInt_t i ; + for (i = 0; i < _theFitter->Result().NTotalParameters(); ++i) { + + TString varName(_theFitter->Result().GetParameterName(i)); + Bool_t isConst(_theFitter->Result().IsParameterFixed(i)) ; + + Double_t xlo = _theFitter->Config().ParSettings(i).LowerLimit(); + Double_t xhi = _theFitter->Config().ParSettings(i).UpperLimit(); + Double_t xerr = _theFitter->Result().Error(i); + Double_t xval = _theFitter->Result().Value(i); + + RooRealVar* var ; + if (varList.getSize()==0) { + + if ((xlosetConstant(isConst) ; + } else { + + var = (RooRealVar*) varList.at(i)->Clone() ; + var->setConstant(isConst) ; + var->setVal(xval) ; + if (xlosetRange(xlo,xhi) ; + } + + if (varName.CompareTo(var->GetName())) { + oocoutI((TObject*)0,Eval) << "RooMinimizerSemiAnalytic::lastMinuitFit: fit parameter '" << varName + << "' stored in variable '" << var->GetName() << "'" << endl ; + } + + } + + if (isConst) { + constPars.addOwned(*var) ; + } else { + var->setError(xerr) ; + floatPars.addOwned(*var) ; + } + } + + res->setConstParList(constPars) ; + res->setInitParList(floatPars) ; + res->setFinalParList(floatPars) ; + res->setMinNLL(_theFitter->Result().MinFcnValue()) ; + res->setEDM(_theFitter->Result().Edm()) ; + res->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ; + res->setStatus(_theFitter->Result().Status()) ; + std::vector globalCC; + TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ; + TMatrixDSym covs(_theFitter->Result().Parameters().size()) ; + for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) { + globalCC.push_back(_theFitter->Result().GlobalCC(ic)); + for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) { + corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii); + covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii); + } + } + res->fillCorrMatrix(globalCC,corrs,covs) ; + + return res; + +} + +#endif diff --git a/src/classes.h b/src/classes.h index 28c454077b2..6acd0d74105 100644 --- a/src/classes.h +++ b/src/classes.h @@ -64,6 +64,9 @@ #include "HiggsAnalysis/CombinedLimit/interface/RooNCSpline_3D_fast.h" #include "HiggsAnalysis/CombinedLimit/interface/RooFuncPdf.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerFcnSemiAnalytic.h" + namespace { struct dictionary { RooBernsteinFast<1> my_RooBernsteinFast_1; diff --git a/src/classes_def.xml b/src/classes_def.xml index 9f72287e6f5..b4d9f36f518 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -208,4 +208,6 @@ + + From aaf434ab60168918b6df72186dab04b0c328e397 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 15 Oct 2019 16:06:21 +0200 Subject: [PATCH 02/16] integration. No derivatives yet. --- interface/CascadeMinimizer.h | 7 +++ src/CascadeMinimizer.cc | 102 +++++++++++++++++++---------------- 2 files changed, 63 insertions(+), 46 deletions(-) diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index 89536f5108c..9644230de30 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -12,6 +12,8 @@ class RooRealVar; #include #include +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" + class CascadeMinimizer { public: @@ -44,6 +46,11 @@ class CascadeMinimizer { private: RooAbsReal & nll_; std::auto_ptr minimizer_; + // + bool isSemiAnalyticMinimizer{false}; + std::auto_ptr minimizerSemiAnalytic_; + std::map derivatives_; + // Mode mode_; static int strategy_; RooRealVar * poi_; diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index fce5753dcb1..19052d4c7ba 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -59,8 +59,18 @@ CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi) void CascadeMinimizer::remakeMinimizer() { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); if (simnll) simnll->setHideRooCategories(true); - minimizer_.reset(); // avoid two copies in memory - minimizer_.reset(new RooMinimizer(nll_)); + + if (runtimedef::get("MINIMIZER_SemiAnalytic")){ + isSemiAnalyticMinimizer=true; + minimizerSemiAnalytic_.reset(); + minimizerSemiAnalytic_.reset(new RooMinimzerSemiAnalytic(nll_,derivatives_)); + } + else{ + isSemiAnalyticMinimizer=false; + minimizer_.reset(); // avoid two copies in memory + minimizer_.reset(new RooMinimizer(nll_)); + } + if (simnll) simnll->setHideRooCategories(false); } @@ -100,29 +110,29 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz simnllbb->setAnalyticBarlowBeeston(true); forceResetMinimizer = true; } - if (forceResetMinimizer || !minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); + if (forceResetMinimizer || !((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); strategy_ = ROOT::Math::MinimizerOptions::DefaultStrategy(); // re-configure - minimizer_->setStrategy(strategy_); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setEps(nominalTol); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); if (approxPreFitTolerance_ > 0) { double tol = std::max(approxPreFitTolerance_, 10. * nominalTol); do { if (verbose > 1) std::cout << "Running pre-fit with " << nominalType << "," << nominalAlgo << " and tolerance " << tol << std::endl; Significance::MinimizerSentry minimizerConfig(nominalType+","+nominalAlgo, tol); - minimizer_->setEps(tol); - minimizer_->setStrategy(approxPreFitStrategy_); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(tol); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(approxPreFitStrategy_); improveOnce(verbose-1, true); if (runtimedef::get("DBG_QUICKEXIT")) { exit(0); } - minimizer_->setEps(nominalTol); - minimizer_->setStrategy(strategy_); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); } while (autoBounds_ && !autoBoundsOk(verbose-1)); } bool outcome; @@ -145,8 +155,8 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz std::cerr << "Will fallback to minimization using " << it->algo << ", strategy " << myStrategy << " and tolerance " << it->tolerance << std::endl; Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Will fallback to minimization using %s, strategy %d and tolerance %g",__LINE__,(it->algo).c_str(),myStrategy,it->tolerance)),Logger::kLogLevelDebug,__func__); } - minimizer_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setStrategy(myStrategy); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(myStrategy); outcome = improveOnce(verbose-2); if (outcome) break; } @@ -171,33 +181,33 @@ bool CascadeMinimizer::improveOnce(int verbose, bool noHesse) bool outcome = false; double tol = ROOT::Math::MinimizerOptions::DefaultTolerance(); static int maxcalls = runtimedef::get("MINIMIZER_MaxCalls"); - if (!minimizer_.get()) remakeMinimizer(); + if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); if (maxcalls) { - minimizer_->setMaxFunctionCalls(maxcalls); - minimizer_->setMaxIterations(maxcalls); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setMaxFunctionCalls(maxcalls); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setMaxIterations(maxcalls); } if (oldFallback_){ - if (optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); - outcome = nllutils::robustMinimize(nll_, *minimizer_, verbose, setZeroPoint_); + if (optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); + if (rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); + outcome = nllutils::robustMinimize(nll_, *((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_), verbose, setZeroPoint_); } else { if (verbose+2>0) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minimisation configured with Type=%s, Algo=%s, strategy=%d, tolerance=%g",__LINE__,myType.c_str(),myAlgo.c_str(),myStrategy,tol)),Logger::kLogLevelInfo,__func__); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if ((!simnll) && optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if ((!simnll) && rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); + if ((!simnll) && optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); + if ((!simnll) && rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); if (firstHesse_ && !noHesse) { - minimizer_->setPrintLevel(std::max(0,verbose-3)); - minimizer_->hesse(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(std::max(0,verbose-3)); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); if (simnll) simnll->updateZeroPoint(); - minimizer_->setPrintLevel(verbose-1); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); } - int status = minimizer_->minimize(myType.c_str(), myAlgo.c_str()); + int status = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(myType.c_str(), myAlgo.c_str()); if (lastHesse_ && !noHesse) { if (simnll) simnll->updateZeroPoint(); - minimizer_->setPrintLevel(std::max(0,verbose-3)); - status = minimizer_->hesse(); - minimizer_->setPrintLevel(verbose-1); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(std::max(0,verbose-3)); + status = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); if (verbose+2>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Hesse finished with status=%d",__LINE__,status)),Logger::kLogLevelDebug,__func__); } if (simnll) simnll->clearZeroPoint(); @@ -231,8 +241,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { utils::setAllConstant(toFreeze, false); remakeMinimizer(); } - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); // for debugging + if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -245,8 +255,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { //TStopwatch tw; // need to re-run Migrad before running minos - minimizer_->minimize(myType.c_str(), "Migrad"); - int iret = minimizer_->minos(params); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(myType.c_str(), "Migrad"); + int iret = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minos(params); if (verbose>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minos finished with status=%d",__LINE__,iret)),Logger::kLogLevelDebug,__func__); //std::cout << "Run Minos in "; tw.Print(); std::cout << std::endl; @@ -270,12 +280,12 @@ bool CascadeMinimizer::hesse(int verbose ) { // Have to reset and minimize again first to get all parameters in remakeMinimizer(); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setEps(nominalTol); - minimizer_->setStrategy(strategy_); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); improveOnce(verbose - 1); } - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); // for debugging + if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -286,7 +296,7 @@ bool CascadeMinimizer::hesse(int verbose ) { } } - int iret = minimizer_->hesse(); + int iret = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); if (setZeroPoint_) { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); @@ -375,9 +385,9 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) bool doMultipleMini = (CascadeMinimizerGlobalConfigs::O().pdfCategories.getSize()>0); if (runtimedef::get(std::string("MINIMIZER_skipDiscreteIterations"))) doMultipleMini=false; // if ( doMultipleMini ) preFit_ = 1; - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-2); - minimizer_->setStrategy(strategy_); + if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-2); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); RooArgSet nuisances = CascadeMinimizerGlobalConfigs::O().nuisanceParameters; @@ -388,13 +398,13 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) freezeDiscParams(true); remakeMinimizer(); - minimizer_->setPrintLevel(verbose-2); - minimizer_->setStrategy(preFit_-1); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-2); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(preFit_-1); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if (optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); - minimizer_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); + if (optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); + if (rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); if (simnll) simnll->clearZeroPoint(); utils::setAllConstant(frozen,false); freezeDiscParams(false); @@ -538,7 +548,7 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, if (hideConstants && simnll) { simnll->setHideConstants(true); if (maskConstraints) simnll->setMaskConstraints(true); - minimizer_.reset(); // will be recreated when needed by whoever needs it + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).reset(); // will be recreated when needed by whoever needs it } std::vector > myCombos; @@ -693,7 +703,7 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, if (hideConstants && simnll) { simnll->setHideConstants(false); if (maskConstraints) simnll->setMaskConstraints(false); - minimizer_.reset(); // will be recreated when needed by whoever needs it + ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).reset(); // will be recreated when needed by whoever needs it } From 767cba7de5d45f62b8cba8b931602dd9ff9a6cc2 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 15 Oct 2019 17:39:25 +0200 Subject: [PATCH 03/16] interface done. initial_fits.sh No derivatives in, yet. --- interface/CascadeMinimizer.h | 3 +- src/CascadeMinimizer.cc | 104 ++++++++++++++++++++--------------- 2 files changed, 60 insertions(+), 47 deletions(-) diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index 9644230de30..fca66012d25 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -14,7 +14,6 @@ class RooRealVar; #include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" - class CascadeMinimizer { public: enum Mode { Constrained, Unconstrained }; @@ -32,7 +31,7 @@ class CascadeMinimizer { RooMinimizer & minimizer() { return *minimizer_; } RooFitResult *save() { return minimizer().save(); } void setStrategy(int strategy) { strategy_ = strategy; } - void setErrorLevel(float errorLevel) { minimizer_->setErrorLevel(errorLevel); } + void setErrorLevel(float errorLevel) { (not isSemiAnalyticMinimizer)? minimizer_->setErrorLevel(errorLevel):minimizerSemiAnalytic_->setErrorLevel(errorLevel); } static void initOptions() ; static void applyOptions(const boost::program_options::variables_map &vm) ; static const boost::program_options::options_description & options() { return options_; } diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index 19052d4c7ba..caed893a7e8 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -15,6 +15,12 @@ #include +namespace local{ + template + class getter{ + }; +}; + boost::program_options::options_description CascadeMinimizer::options_("Cascade Minimizer options"); std::vector CascadeMinimizer::fallbacks_; bool CascadeMinimizer::preScan_; @@ -63,7 +69,7 @@ void CascadeMinimizer::remakeMinimizer() { if (runtimedef::get("MINIMIZER_SemiAnalytic")){ isSemiAnalyticMinimizer=true; minimizerSemiAnalytic_.reset(); - minimizerSemiAnalytic_.reset(new RooMinimzerSemiAnalytic(nll_,derivatives_)); + minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_,derivatives_)); } else{ isSemiAnalyticMinimizer=false; @@ -74,6 +80,7 @@ void CascadeMinimizer::remakeMinimizer() { if (simnll) simnll->setHideRooCategories(false); } + bool CascadeMinimizer::freezeDiscParams(const bool freeze) { static bool freezeDisassParams = runtimedef::get(std::string("MINIMIZER_freezeDisassociatedParams")); @@ -110,29 +117,29 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz simnllbb->setAnalyticBarlowBeeston(true); forceResetMinimizer = true; } - if (forceResetMinimizer || !((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); + if (forceResetMinimizer || !(minimizer_.get() or minimizerSemiAnalytic_.get())) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); strategy_ = ROOT::Math::MinimizerOptions::DefaultStrategy(); // re-configure - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); if (approxPreFitTolerance_ > 0) { double tol = std::max(approxPreFitTolerance_, 10. * nominalTol); do { if (verbose > 1) std::cout << "Running pre-fit with " << nominalType << "," << nominalAlgo << " and tolerance " << tol << std::endl; Significance::MinimizerSentry minimizerConfig(nominalType+","+nominalAlgo, tol); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(tol); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(approxPreFitStrategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(approxPreFitStrategy_):minimizerSemiAnalytic_->setStrategy(approxPreFitStrategy_); improveOnce(verbose-1, true); if (runtimedef::get("DBG_QUICKEXIT")) { exit(0); } - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(approxPreFitStrategy_):minimizerSemiAnalytic_->setStrategy(strategy_); } while (autoBounds_ && !autoBoundsOk(verbose-1)); } bool outcome; @@ -155,8 +162,8 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz std::cerr << "Will fallback to minimization using " << it->algo << ", strategy " << myStrategy << " and tolerance " << it->tolerance << std::endl; Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Will fallback to minimization using %s, strategy %d and tolerance %g",__LINE__,(it->algo).c_str(),myStrategy,it->tolerance)),Logger::kLogLevelDebug,__func__); } - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(myStrategy); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()):minimizerSemiAnalytic_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(myStrategy):minimizerSemiAnalytic_->setStrategy(myStrategy); outcome = improveOnce(verbose-2); if (outcome) break; } @@ -181,33 +188,38 @@ bool CascadeMinimizer::improveOnce(int verbose, bool noHesse) bool outcome = false; double tol = ROOT::Math::MinimizerOptions::DefaultTolerance(); static int maxcalls = runtimedef::get("MINIMIZER_MaxCalls"); - if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); if (maxcalls) { - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setMaxFunctionCalls(maxcalls); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setMaxIterations(maxcalls); + (not isSemiAnalyticMinimizer)?minimizer_->setMaxFunctionCalls(maxcalls):minimizerSemiAnalytic_->setMaxFunctionCalls(maxcalls); + (not isSemiAnalyticMinimizer)?minimizer_->setMaxIterations(maxcalls):minimizerSemiAnalytic_->setMaxIterations(maxcalls); } if (oldFallback_){ - if (optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); - outcome = nllutils::robustMinimize(nll_, *((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_), verbose, setZeroPoint_); + if (isSemiAnalyticMinimizer)throw std::logic_error("unimplemented");// check robust minimizer stuff with Minimizer + + if (optConst){ + minimizer_->optimizeConst(std::max(0,optConst)); + } + if (rooFitOffset) { minimizer_->setOffsetting(std::max(0,rooFitOffset));} + + outcome = nllutils::robustMinimize(nll_, *minimizer_, verbose, setZeroPoint_); } else { if (verbose+2>0) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minimisation configured with Type=%s, Algo=%s, strategy=%d, tolerance=%g",__LINE__,myType.c_str(),myAlgo.c_str(),myStrategy,tol)),Logger::kLogLevelInfo,__func__); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if ((!simnll) && optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); - if ((!simnll) && rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); + if ((!simnll) && optConst) (not isSemiAnalyticMinimizer)?minimizer_->optimizeConst(std::max(0,optConst)):minimizerSemiAnalytic_->optimizeConst(std::max(0,optConst)); + if ((!simnll) && rooFitOffset) (not isSemiAnalyticMinimizer)?minimizer_->setOffsetting(std::max(0,rooFitOffset)):minimizerSemiAnalytic_->setOffsetting(std::max(0,rooFitOffset)); if (firstHesse_ && !noHesse) { - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(std::max(0,verbose-3)); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(std::max(0,verbose-3)):minimizerSemiAnalytic_->setPrintLevel(std::max(0,verbose-3)); + (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); if (simnll) simnll->updateZeroPoint(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); } - int status = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(myType.c_str(), myAlgo.c_str()); + int status = (not isSemiAnalyticMinimizer) ? minimizer_->minimize(myType.c_str(), myAlgo.c_str()):minimizerSemiAnalytic_->minimize(myType.c_str(), myAlgo.c_str()); if (lastHesse_ && !noHesse) { if (simnll) simnll->updateZeroPoint(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(std::max(0,verbose-3)); - status = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(std::max(0,verbose-3)):minimizerSemiAnalytic_->setPrintLevel(std::max(0,verbose-3)); + status = (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); if (verbose+2>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Hesse finished with status=%d",__LINE__,status)),Logger::kLogLevelDebug,__func__); } if (simnll) simnll->clearZeroPoint(); @@ -241,8 +253,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { utils::setAllConstant(toFreeze, false); remakeMinimizer(); } - if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); // for debugging + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -255,8 +267,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { //TStopwatch tw; // need to re-run Migrad before running minos - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(myType.c_str(), "Migrad"); - int iret = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minos(params); + (not isSemiAnalyticMinimizer)?minimizer_->minimize(myType.c_str(), "Migrad"):minimizerSemiAnalytic_->minimize(myType.c_str(), "Migrad"); + int iret = (not isSemiAnalyticMinimizer)?minimizer_->minos(params):minimizerSemiAnalytic_->minos(params); if (verbose>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minos finished with status=%d",__LINE__,iret)),Logger::kLogLevelDebug,__func__); //std::cout << "Run Minos in "; tw.Print(); std::cout << std::endl; @@ -280,12 +292,12 @@ bool CascadeMinimizer::hesse(int verbose ) { // Have to reset and minimize again first to get all parameters in remakeMinimizer(); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setEps(nominalTol); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); improveOnce(verbose - 1); } - if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-1); // for debugging + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -296,7 +308,7 @@ bool CascadeMinimizer::hesse(int verbose ) { } } - int iret = ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->hesse(); + int iret = (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); if (setZeroPoint_) { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); @@ -385,9 +397,9 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) bool doMultipleMini = (CascadeMinimizerGlobalConfigs::O().pdfCategories.getSize()>0); if (runtimedef::get(std::string("MINIMIZER_skipDiscreteIterations"))) doMultipleMini=false; // if ( doMultipleMini ) preFit_ = 1; - if (!((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).get()) remakeMinimizer(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-2); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(strategy_); + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-2):minimizerSemiAnalytic_->setPrintLevel(verbose-2); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); RooArgSet nuisances = CascadeMinimizerGlobalConfigs::O().nuisanceParameters; @@ -398,13 +410,13 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) freezeDiscParams(true); remakeMinimizer(); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setPrintLevel(verbose-2); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setStrategy(preFit_-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-2):minimizerSemiAnalytic_->setPrintLevel(verbose-2); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(preFit_-1):minimizerSemiAnalytic_->setStrategy(preFit_-1); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if (optConst) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->setOffsetting(std::max(0,rooFitOffset)); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_)->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); + if (optConst) (not isSemiAnalyticMinimizer)?minimizer_->optimizeConst(std::max(0,optConst)):minimizerSemiAnalytic_->optimizeConst(std::max(0,optConst)); + if (rooFitOffset) (not isSemiAnalyticMinimizer)?minimizer_->setOffsetting(std::max(0,rooFitOffset)):minimizerSemiAnalytic_->setOffsetting(std::max(0,rooFitOffset)); + (not isSemiAnalyticMinimizer)?minimizer_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()):minimizerSemiAnalytic_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); if (simnll) simnll->clearZeroPoint(); utils::setAllConstant(frozen,false); freezeDiscParams(false); @@ -548,7 +560,8 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, if (hideConstants && simnll) { simnll->setHideConstants(true); if (maskConstraints) simnll->setMaskConstraints(true); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).reset(); // will be recreated when needed by whoever needs it + minimizer_.reset(); + minimizerSemiAnalytic_.reset(); // will be recreated when needed by whoever needs it } std::vector > myCombos; @@ -703,7 +716,8 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, if (hideConstants && simnll) { simnll->setHideConstants(false); if (maskConstraints) simnll->setMaskConstraints(false); - ((isSemiAnalyticMinimizer)?minimizerSemiAnalytic_:minimizer_).reset(); // will be recreated when needed by whoever needs it + minimizer_.reset(); // will be recreated when needed by whoever needs it + minimizerSemiAnalytic_.reset(); } From 863dfd8da476d6980fcb4e3599dc9e9cdc9a5221 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 15 Oct 2019 19:32:59 +0200 Subject: [PATCH 04/16] small changes --- interface/CascadeMinimizer.h | 3 ++- interface/RooMinimizerFcnSemiAnalytic.h | 2 +- src/CascadeMinimizer.cc | 6 ++++-- src/RooMinimizerFcnSemiAnalytic.cxx | 6 +++--- src/RooMinimizerSemiAnalytic.cc | 10 +++++----- 5 files changed, 15 insertions(+), 12 deletions(-) diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index fca66012d25..5e67b46f528 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -13,6 +13,7 @@ class RooRealVar; #include #include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include class CascadeMinimizer { public: @@ -28,7 +29,7 @@ class CascadeMinimizer { bool improve(int verbose=0, bool cascade=true, bool forceResetMinimizer=false); // declare nuisance parameters for pre-fit void setNuisanceParameters(const RooArgSet *nuis) { nuisances_ = nuis; } - RooMinimizer & minimizer() { return *minimizer_; } + RooMinimizer & minimizer() { if(isSemiAnalyticMinimizer) throw std::runtime_error("unimplemented"); return *minimizer_; } RooFitResult *save() { return minimizer().save(); } void setStrategy(int strategy) { strategy_ = strategy; } void setErrorLevel(float errorLevel) { (not isSemiAnalyticMinimizer)? minimizer_->setErrorLevel(errorLevel):minimizerSemiAnalytic_->setErrorLevel(errorLevel); } diff --git a/interface/RooMinimizerFcnSemiAnalytic.h b/interface/RooMinimizerFcnSemiAnalytic.h index a885fd27cfc..dac08c85e51 100644 --- a/interface/RooMinimizerFcnSemiAnalytic.h +++ b/interface/RooMinimizerFcnSemiAnalytic.h @@ -40,7 +40,7 @@ class RooMinimizerFcnSemiAnalytic : RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other); virtual ~RooMinimizerFcnSemiAnalytic(); - ROOT::Math::IBaseFunctionMultiDim* Clone() const override; + ROOT::Math::IMultiGradFunction* Clone() const override; unsigned int NDim() const override{ return _nDim; } RooArgList* GetFloatParamList() { return _floatParamList; } diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index caed893a7e8..cac4c4078c5 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -66,14 +66,16 @@ void CascadeMinimizer::remakeMinimizer() { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); if (simnll) simnll->setHideRooCategories(true); + minimizer_.reset(); // avoid two copies in memory + minimizerSemiAnalytic_.reset(); + if (runtimedef::get("MINIMIZER_SemiAnalytic")){ isSemiAnalyticMinimizer=true; - minimizerSemiAnalytic_.reset(); minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_,derivatives_)); + //minimizerSemiAnalytic_->setVerbose(kTRUE); // DEBUG } else{ isSemiAnalyticMinimizer=false; - minimizer_.reset(); // avoid two copies in memory minimizer_.reset(new RooMinimizer(nll_)); } diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx index 62e9082598c..d4d01d1f67f 100644 --- a/src/RooMinimizerFcnSemiAnalytic.cxx +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -125,7 +125,7 @@ RooMinimizerFcnSemiAnalytic::~RooMinimizerFcnSemiAnalytic() } -ROOT::Math::IBaseFunctionMultiDim* RooMinimizerFcnSemiAnalytic::Clone() const +ROOT::Math::IMultiGradFunction* RooMinimizerFcnSemiAnalytic::Clone() const { return new RooMinimizerFcnSemiAnalytic(*this) ; } @@ -563,7 +563,7 @@ double RooMinimizerFcnSemiAnalytic::DoEval(const double *x) const Bool_t first(kTRUE) ; /*ooccoutW(_context,Minimization)*/std::cout << "Parameter values: " ; - //for (const auto par : *_floatParamList) { //new loop + //for (const auto par : *_floatParamList) { //new loop} for (int ip =0 ;ip < _floatParamList->getSize();++ip){ const auto par=_floatParamList->at(ip); auto var = static_cast(par); @@ -644,7 +644,7 @@ double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int i ROOT use this has step above, maybe copy it dynamically? TODO double MultiNumGradFunction::fgEps = 0.001; - double MultiNumGradFunction::DoDerivative (const double * x, unsigned int icoord ) const { + double MultiNumGradFunction::DoDerivative (const double * x, unsigned int icoord ) const // calculate derivative using mathcore derivator class // step size can be changes using SetDerivPrecision() diff --git a/src/RooMinimizerSemiAnalytic.cc b/src/RooMinimizerSemiAnalytic.cc index 08f2dae6224..9b14684e24f 100644 --- a/src/RooMinimizerSemiAnalytic.cc +++ b/src/RooMinimizerSemiAnalytic.cc @@ -309,7 +309,7 @@ Int_t RooMinimizerSemiAnalytic::minimize(const char* type, const char* alg) RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; RooAbsReal::clearEvalErrorLog() ; - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -338,7 +338,7 @@ Int_t RooMinimizerSemiAnalytic::migrad() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -501,7 +501,7 @@ Int_t RooMinimizerSemiAnalytic::seek() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -530,7 +530,7 @@ Int_t RooMinimizerSemiAnalytic::simplex() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -559,7 +559,7 @@ Int_t RooMinimizerSemiAnalytic::improve() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; From 374497fc46bf3afeee7b77a56e18a833d386de37 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 15 Oct 2019 19:37:54 +0200 Subject: [PATCH 05/16] debug info --- src/RooMinimizerFcnSemiAnalytic.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx index d4d01d1f67f..8648e9199fd 100644 --- a/src/RooMinimizerFcnSemiAnalytic.cxx +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -84,7 +84,9 @@ RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooM _nDim = _floatParamList->getSize(); + _verbose=true; // print the variables for which we have derivatives? updateFloatVec() ; + _verbose=verbose; // Save snapshot of initial lists _initFloatParamList = (RooArgList*) _floatParamList->snapshot(kFALSE) ; @@ -535,6 +537,7 @@ void RooMinimizerFcnSemiAnalytic::updateFloatVec() double RooMinimizerFcnSemiAnalytic::DoEval(const double *x) const { + //std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "DoEval"< Date: Thu, 17 Sep 2020 16:41:06 +0200 Subject: [PATCH 06/16] last things --- interface/CachingNLL.h | 5 +++++ interface/ProcessNormalization.h | 3 +++ interface/RooMinimizerFcnSemiAnalytic.h | 1 + interface/RooMinimizerSemiAnalytic.h | 1 + interface/SimpleGaussianConstraint.h | 2 ++ src/CachingNLL.cc | 2 +- 6 files changed, 13 insertions(+), 1 deletion(-) diff --git a/interface/CachingNLL.h b/interface/CachingNLL.h index 35578f4fa31..4d7f2dd25eb 100644 --- a/interface/CachingNLL.h +++ b/interface/CachingNLL.h @@ -21,6 +21,7 @@ #include class RooMultiPdf; +class SimNLLDerivativesHelper; // Part zero: ArgSet checker namespace cacheutils { @@ -112,6 +113,7 @@ class OptimizedCachingPdfT : public CachingPdf { CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; class CachingAddNLL : public RooAbsReal { + friend SimNLLDerivativesHelper; // probably not needed w/ data public: CachingAddNLL(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data, bool includeZeroWeights = false) ; CachingAddNLL(const CachingAddNLL &other, const char *name = 0) ; @@ -125,6 +127,7 @@ class CachingAddNLL : public RooAbsReal { virtual RooArgSet* getParameters(const RooArgSet* depList, Bool_t stripDisconnected = kTRUE) const ; double sumWeights() const { return sumWeights_; } const RooAbsPdf *pdf() const { return pdf_; } + const RooAbsData *data() const {return data_;} void setZeroPoint() ; void clearZeroPoint() ; void clearConstantZeroPoint() ; @@ -186,6 +189,8 @@ class CachingSimNLL : public RooAbsReal { // trap this call, since we don't care about propagating it to the sub-components virtual void constOptimizeTestStatistic(ConstOpCode opcode, Bool_t doAlsoTrackingOpt=kTRUE) { } private: + friend SimNLLDerivativesHelper; + void setup_(); RooSimultaneous *pdfOriginal_; const RooAbsData *dataOriginal_; diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 0721e198f25..91f992cb19c 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -14,6 +14,8 @@ ProcessNormalization is helper class for implementing process normalizations END_HTML */ // +class SimNLLDerivativesHelper; + class ProcessNormalization : public RooAbsReal { public: ProcessNormalization() : nominalValue_(1) {} @@ -33,6 +35,7 @@ class ProcessNormalization : public RooAbsReal { Double_t evaluate() const; private: + friend SimNLLDerivativesHelper; // ---- PERSISTENT ---- double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas diff --git a/interface/RooMinimizerFcnSemiAnalytic.h b/interface/RooMinimizerFcnSemiAnalytic.h index dac08c85e51..16377ea3783 100644 --- a/interface/RooMinimizerFcnSemiAnalytic.h +++ b/interface/RooMinimizerFcnSemiAnalytic.h @@ -2,6 +2,7 @@ * Authors: * * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * * Code based on the equivalent roofit code/file. * + * Wed Oct 16 14:04:34 CEST 2019 * * * Redistribution and use in source and binary forms, * * with or without modification, are permitted according to the terms * diff --git a/interface/RooMinimizerSemiAnalytic.h b/interface/RooMinimizerSemiAnalytic.h index 2ed1857d3c5..063d86447f4 100644 --- a/interface/RooMinimizerSemiAnalytic.h +++ b/interface/RooMinimizerSemiAnalytic.h @@ -2,6 +2,7 @@ * Authors: * * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * * Code based on the equivalent roofit code/file. * + * Wed Oct 16 14:04:34 CEST 2019 * * * * * Redistribution and use in source and binary forms, * diff --git a/interface/SimpleGaussianConstraint.h b/interface/SimpleGaussianConstraint.h index f1f2fd4c2dc..364dd271c5a 100644 --- a/interface/SimpleGaussianConstraint.h +++ b/interface/SimpleGaussianConstraint.h @@ -17,6 +17,8 @@ class SimpleGaussianConstraint : public RooGaussian { inline virtual ~SimpleGaussianConstraint() { } const RooAbsReal & getX() const { return x.arg(); } + const RooAbsReal & getMean() const { return mean.arg(); } + const RooAbsReal & getSigma() const { return sigma.arg(); } double getLogValFast() const { if (_valueDirty) { diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index 23952be84cd..93c1b90ecdd 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -611,7 +611,7 @@ cacheutils::CachingAddNLL::evaluate() const boost::ptr_vector::iterator itp = pdfs_.begin();//, edp = pdfs_.end(); std::vector::const_iterator itw; //bgw = weights_.begin();//, edw = weights_.end(); std::vector::iterator its, bgs = partialSum_.begin(), eds = partialSum_.end(); - double sumCoeff = 0; + double sumCoeff = 0.; bool allBasicIntegralsOk = (basicIntegrals_ == 1); //std::cout << "Performing evaluation of " << GetName() << std::endl; for ( ; itc != edc; ++itp, ++itc ) { From b61e1e8369e59e60d6b2aa70535263833f6bd897 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Thu, 1 Jul 2021 19:02:58 +0200 Subject: [PATCH 07/16] first version of derivatives in semianalytical gradient --- interface/CachingNLL.h | 2 + interface/CascadeMinimizer.h | 2 +- interface/ProcessNormalization.h | 2 + interface/SimNLLDerivativesHelper.h | 61 ++++++++++ src/CachingNLL.cc | 1 - src/CascadeMinimizer.cc | 10 +- src/SimNLLDerivativesHelper.cc | 180 ++++++++++++++++++++++++++++ 7 files changed, 255 insertions(+), 3 deletions(-) create mode 100644 interface/SimNLLDerivativesHelper.h create mode 100644 src/SimNLLDerivativesHelper.cc diff --git a/interface/CachingNLL.h b/interface/CachingNLL.h index 4d7f2dd25eb..84e2806a0ab 100644 --- a/interface/CachingNLL.h +++ b/interface/CachingNLL.h @@ -22,6 +22,7 @@ class RooMultiPdf; class SimNLLDerivativesHelper; +class DerivativeLogNormal; // Part zero: ArgSet checker namespace cacheutils { @@ -114,6 +115,7 @@ CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; class CachingAddNLL : public RooAbsReal { friend SimNLLDerivativesHelper; // probably not needed w/ data + friend DerivativeLogNormal; public: CachingAddNLL(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data, bool includeZeroWeights = false) ; CachingAddNLL(const CachingAddNLL &other, const char *name = 0) ; diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index 5e67b46f528..e6575c7fd8e 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -49,7 +49,7 @@ class CascadeMinimizer { // bool isSemiAnalyticMinimizer{false}; std::auto_ptr minimizerSemiAnalytic_; - std::map derivatives_; + std::map derivatives_; // not used // Mode mode_; static int strategy_; diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 91f992cb19c..61bf610f216 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -15,6 +15,7 @@ END_HTML */ // class SimNLLDerivativesHelper; +class DerivativeLogNormal; class ProcessNormalization : public RooAbsReal { public: @@ -36,6 +37,7 @@ class ProcessNormalization : public RooAbsReal { private: friend SimNLLDerivativesHelper; + friend DerivativeLogNormal; // ---- PERSISTENT ---- double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h new file mode 100644 index 00000000000..7a44cced7b3 --- /dev/null +++ b/interface/SimNLLDerivativesHelper.h @@ -0,0 +1,61 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, CERN (CH), andrea.carlo.marini@cern.ch * + * Code based on roofit code. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef SimNLLDerivHelper_H +#define SimNLLDerivHelper_H +#include +#include +#include +#include +#include "RooAddition.h" +#include "RooDataSet.h" +#include "HiggsAnalysis/CombinedLimit/interface/CachingNLL.h" +#include "HiggsAnalysis/CombinedLimit/interface/ProcessNormalization.h" + +/* This class takes in input the CachingSimNLL and returns the derivatives of what it knows + * + */ + +class DerivativeLogNormal: public RooAbsReal +{ + const RooDataSet * data_{nullptr}; // not owned. + const cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL + std::string kappaname_{""}; + + public: + DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& kappaname); + ~DerivativeLogNormal(){}; + virtual Bool_t isDerived() const { return kTRUE; } + virtual DerivativeLogNormal *clone(const char *name = 0) const ; + const RooDataSet *data() const {return data_;} + const cacheutils::CachingAddNLL *pdf() const { return pdf_; } + + virtual Double_t evaluate() const ; // cacheutils::ReminderSum + + // name-> index association per process + std::vector kappa_pos_; +}; + + +class SimNLLDerivativesHelper +{ + public: + SimNLLDerivativesHelper( cacheutils::CachingSimNLL * nll ){nll_=nll;} + ~SimNLLDerivativesHelper(){}; + void init(); + bool verbose{true}; // debug + std::map derivatives_; + //std::map& derivatives(){return derivatives;} + private: + cacheutils::CachingSimNLL* nll_{nullptr}; //not owned, keep pointer +}; + +#endif + diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index 93c1b90ecdd..c8194eca339 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -985,7 +985,6 @@ cacheutils::CachingSimNLL::setup_() factorizedPdf_.release(); simpdf = dynamic_cast(pdfclone); } - std::auto_ptr catClone((RooAbsCategoryLValue*) simpdf->indexCat().Clone()); pdfs_.resize(catClone->numBins(NULL), 0); diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index 6bec59f86c9..6c887fb1f8d 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -15,6 +15,8 @@ #include +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" + namespace local{ template class getter{ @@ -71,7 +73,13 @@ void CascadeMinimizer::remakeMinimizer() { if (runtimedef::get("MINIMIZER_SemiAnalytic")){ isSemiAnalyticMinimizer=true; - minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_,derivatives_)); + if (simnll==nullptr) { + std::cout<<"[ERROR]::[CascadeMinimizer] I need simnll to have SemiAnalytic derivatives done by SimNLLDerivativesHelper"<setVerbose(kTRUE); // DEBUG } else{ diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc new file mode 100644 index 00000000000..7ded39bf956 --- /dev/null +++ b/src/SimNLLDerivativesHelper.cc @@ -0,0 +1,180 @@ +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" + +DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& kappaname) : + RooAbsReal(name, title), + data_(data), + pdf_(pdf), + kappaname_(kappaname) +{ + for (unsigned int i=0;ipdfs_.size();++i) + { + int val=-1; + ProcessNormalization *p = dynamic_cast (pdf_->coeffs_[i]); + for (unsigned int j=0 ; jlogKappa_.size();++j) + { + if ( p->thetaList_.at(j)->GetName() == kappaname_ ) {val=int(j);} + } + if (val==-1) { + std::cout<<"[ERROR]: Unable to find kappa corresponding to "<GetName() < sum of expectations in the bin + * lambdat_b -> sum of expecations times logK + */ + + double sum=0.0; + // caching: TODO + for (int ib=0;ibnumEntries();++ib) + { + data_->get(ib) ; + double db = data_->weight(); + double lambda=pdf_->getVal(); + double lambdat=0.0; + for (unsigned int i=0;ipdfs_.size();++i) // loop over processes + { + ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); + double logK = c -> logKappa_ [kappa_pos_[i]]; + lambdat+= c->getVal() * pdf_->pdfs_.at(i).pdf()->getVal() * logK; + } + sum += lambdat - db*lambdat/lambda; + } + return sum; + +} + +// look inside and check the parameters associated to kappas +// compute the sums +// construct RooAbsReal +// keep track of what I can't do like that. + +void SimNLLDerivativesHelper::init(){ + if (verbose) std::cout<<"Init SimNLLDerivativesHelper"<Delete(); + derivatives_.clear(); + //--- + + std::set logNormal; + // Find LogNormal candidates + for (auto f : nll_->constrainPdfsFast_) + { + if (verbose) std::cout<<"inserting candidate: "<< f->getX().GetName()<getX().GetName() ); + } + + // loop over sim components + unsigned idx=0; + for (std::vector::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + bool isWeighted = data->isWeighted(); // binned vs unbinned?!? TBC + + for(unsigned proc =0 ; proc < pdf ->pdfs_.size();++proc) + { + RooAbsReal * coeff = pdf -> coeffs_[proc] ; + ProcessNormalization * pn = dynamic_cast (coeff); + if (pn == nullptr or not isWeighted) { // remove all the lognormal candidates. Don't know how to deal with them + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing from the list of known derivatives because" + << ((pn==nullptr)? std::string(" coeff is not processNormalization but")+ std::string( coeff->ClassName() ): "") + << ((pn==nullptr and not isWeighted) ? " and":"" ) + << ((not isWeighted)?" dataset is not weighted (unbinned?)":"") + << "the following variables" + <getVariables()); + for (int idx=0; idxGetName()) != logNormal.end()) { + logNormal.erase(v->GetName()); + if (verbose) std::cout<<"|"<GetName(); + } + } + //-- + l = RooArgList(*pdf->getVariables()); + for (int idx=0; idxGetName()) != logNormal.end()) { + logNormal.erase(v->GetName()); + if (verbose) std::cout<<"|"<GetName(); + } + } + continue; + } // not a process normalization coefficient + + // remove all asymmThetaList + //for (auto v : pn->asymmThetaList_) + for(int idx=0; idx< pn->asymmThetaList_.getSize() ;++idx){ + RooAbsArg*v=pn->asymmThetaList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<<" from the list of known derivatives because appears in an asymmThetaList " + <GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + } + // remove all otherFactorList + //for (auto v : pn->otherFactorList_) + for(int idx=0; idx< pn->otherFactorList_.getSize() ;++idx){ + RooAbsArg*v=pn->otherFactorList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<<" from the list of known derivatives because appears in an otherFactors " + <GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + } + + } + + // cross check derivatives with logNormal + //for (auto d : derivatives){ + // if (logNormal.find(d.first) == logNormal.end()){ + // //remove derivative + // derivatives.erase(d.first()); + // } + //} + + } + + for (auto name : logNormal){ + RooArgList list; + // loop over sim components + unsigned idx=0; + for (std::vector::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + // construct the derivative as sum + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string kappaname) + DerivativeLogNormal der("","",pdf,dynamic_cast(data),name); + list.add(der); // or addOwned? + } + // add derivative of constraint + // nll_->constrainPdfsFast_ + for(unsigned i=0;iconstrainPdfsFast_.size();++i) + { + if (nll_->constrainPdfsFast_[i]->getX().GetName() != name) continue; + RooFormulaVar constr("constr","(@0-@1)/(@2*@2)",RooArgList(nll_->constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma()));// (x-mean)/sigma^2 + list.add(constr) ; // or addOwned + } + derivatives_[name] = new RooAddition(("derivative_"+name).c_str(),"",list); + } + + +} + From ced2ea7011e683a0692ecf5828e4517d63477faa Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Fri, 2 Jul 2021 17:09:54 +0200 Subject: [PATCH 08/16] work in progress --- interface/SimNLLDerivativesHelper.h | 4 +- src/RooMinimizerFcnSemiAnalytic.cxx | 7 ++ src/SimNLLDerivativesHelper.cc | 140 ++++++++++++++++++++++++---- 3 files changed, 130 insertions(+), 21 deletions(-) diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h index 7a44cced7b3..ee0814a1231 100644 --- a/interface/SimNLLDerivativesHelper.h +++ b/interface/SimNLLDerivativesHelper.h @@ -30,7 +30,7 @@ class DerivativeLogNormal: public RooAbsReal std::string kappaname_{""}; public: - DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& kappaname); + DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& kappaname,int&found); ~DerivativeLogNormal(){}; virtual Bool_t isDerived() const { return kTRUE; } virtual DerivativeLogNormal *clone(const char *name = 0) const ; @@ -41,6 +41,8 @@ class DerivativeLogNormal: public RooAbsReal // name-> index association per process std::vector kappa_pos_; + + bool verbose{true}; }; diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx index 8648e9199fd..12277dc2815 100644 --- a/src/RooMinimizerFcnSemiAnalytic.cxx +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -618,6 +618,13 @@ double RooMinimizerFcnSemiAnalytic::DoDerivative(const double * x, unsigned int throw std::runtime_error(Form("Derivative vector is too small. (requested coordinate) %u < (derivative size) %lu",icoord,_derivParamVec.size())); } if (_derivParamVec[icoord] == nullptr) return DoNumericalDerivative(x,icoord); + + //verify correctness of analytical derivatives // DEBUG + if (true){ + std::cout<<"[RooMinimizerFcnSemiAnalytic][DEBUG]"<< "Doing derivative for "<getVal()<getVal(); diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index 7ded39bf956..03bf61d7e7b 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -1,33 +1,58 @@ #include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" -DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& kappaname) : +DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& kappaname,int&found) : RooAbsReal(name, title), data_(data), pdf_(pdf), kappaname_(kappaname) { - for (unsigned int i=0;ipdfs_.size();++i) + found=0; // keep track if at least one depends on the given kappa + for (unsigned int i=0;ipdfs_.size();++i) // process loop { + if(verbose) std::cout<<"[DerivativeLogNormal]: loop "<< i <<"/"<pdfs_.size()<<" considering pdf "<GetName() < (pdf_->coeffs_[i]); - for (unsigned int j=0 ; jlogKappa_.size();++j) - { - if ( p->thetaList_.at(j)->GetName() == kappaname_ ) {val=int(j);} + if (verbose) std::cout<<"[DerivativeLogNormal]: coeff is of class "<coeffs_[i]->ClassName()<(pdf_->coeffs_[i]); + if (p == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) p = dynamic_cast( &pp->components()[ip]); + } + } + // + if (p != nullptr) { // unable to add the sum for this process + + for (unsigned int j=0 ; jlogKappa_.size();++j) + { + if ( p->thetaList_.at(j)->GetName() == kappaname_ ) {val=int(j);} + } + if (val==-1) { + if(verbose) std::cout<<"[DerivativeLogNormal]: Unable to find kappa corresponding to "<GetName() <GetName() << "because not ProcessNormalization" <GetName() <GetName() < sum of expectations in the bin @@ -45,8 +70,18 @@ Double_t DerivativeLogNormal::evaluate() const { for (unsigned int i=0;ipdfs_.size();++i) // loop over processes { ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); - double logK = c -> logKappa_ [kappa_pos_[i]]; - lambdat+= c->getVal() * pdf_->pdfs_.at(i).pdf()->getVal() * logK; + + // fix RooProduct + RooProduct *pp = dynamic_cast(pdf_->coeffs_[i]); + if (c == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); + } + } + ///--- + double logK = (kappa_pos_[i]>=0)? c -> logKappa_ [kappa_pos_[i]] : 0.0; + lambdat+= pdf_->coeffs_[i]->getVal() * pdf_->pdfs_.at(i).pdf()->getVal() * logK; } sum += lambdat - db*lambdat/lambda; } @@ -74,7 +109,14 @@ void SimNLLDerivativesHelper::init(){ logNormal.insert( f->getX().GetName() ); } + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][init] List of constraint pdfs parameters:"<::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { cacheutils::CachingAddNLL * pdf = *it; @@ -82,20 +124,58 @@ void SimNLLDerivativesHelper::init(){ const RooAbsData* data= pdf->data(); bool isWeighted = data->isWeighted(); // binned vs unbinned?!? TBC + if (verbose) { std::cout << "[SimNLLDerivativesHelper][init] > considering pdf "<GetName()<< ((isWeighted)?" isWeighted":" notWeighted") < toRemove; + for(unsigned proc =0 ; proc < pdf ->pdfs_.size();++proc) { - RooAbsReal * coeff = pdf -> coeffs_[proc] ; + RooAbsReal * coeff = pdf -> coeffs_[proc] ; // what to do for RooProduct? + ProcessNormalization * pn = dynamic_cast (coeff); + RooProduct *pp = dynamic_cast(coeff); + + if (pn == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) pn = dynamic_cast( &pp->components()[ip]); + else { + if (verbose) { std::cout <<"components:";} + RooArgList l = RooArgList(*pp->components()[ip].getVariables()); + for (int idx=0; idxGetName()) != logNormal.end()) { + logNormal.erase(v->GetName()); + if (verbose) std::cout<<"|"<GetName(); + } + } + if (verbose) { std::cout <ClassName() ): "") + << ((pn==nullptr)? std::string(" coeff is not processNormalization but ")+ std::string(coeff->ClassName()): "") << ((pn==nullptr and not isWeighted) ? " and":"" ) << ((not isWeighted)?" dataset is not weighted (unbinned?)":"") - << "the following variables" <Print("V"); + std::cout<<"-----------------"<getVariables()); for (int idx=0; idxgetVariables()); for (int idx=0; idxGetName(); } } + if (verbose) { std::cout<> Removing all asymmThetaList "<asymmThetaList_) for(int idx=0; idx< pn->asymmThetaList_.getSize() ;++idx){ @@ -129,6 +212,7 @@ void SimNLLDerivativesHelper::init(){ } // remove all otherFactorList //for (auto v : pn->otherFactorList_) + std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all otherfactorsList "<otherFactorList_.getSize() ;++idx){ RooAbsArg*v=pn->otherFactorList_.at(idx); if (verbose) { @@ -139,6 +223,8 @@ void SimNLLDerivativesHelper::init(){ if (logNormal.find(v->GetName()) != logNormal.end()) logNormal.erase(v->GetName()); } + std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Continue loop "< Continue loop "<data(); //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string kappaname) - DerivativeLogNormal der("","",pdf,dynamic_cast(data),name); - list.add(der); // or addOwned? + int found; + DerivativeLogNormal *der= new DerivativeLogNormal("","",pdf,dynamic_cast(data),name,found); + if(found) list.addOwned(*der); // or addOwned or add? + else der->Delete(); } // add derivative of constraint // nll_->constrainPdfsFast_ + std::cout<<"[SimNLLDerivativesHelper][init]"<< "Adding constraint term for "<< name <constrainPdfsFast_.size();++i) { + if(verbose) std::cout<<"[SimNLLDerivativesHelper][init]"<< " Considering "<< nll_->constrainPdfsFast_[i]->getX().GetName() <constrainPdfsFast_[i]->getX().GetName() != name) continue; - RooFormulaVar constr("constr","(@0-@1)/(@2*@2)",RooArgList(nll_->constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma()));// (x-mean)/sigma^2 - list.add(constr) ; // or addOwned + + if(verbose) std::cout<<"[SimNLLDerivativesHelper][init]"<< "Loop. found constraint "<< nll_->constrainPdfsFast_[i]->getX().GetName() <<" == "<constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma()));// (x-mean)/sigma^2 + list.addOwned(*constr) ; // or addOwned or add? + } + if (verbose) {std::cout <<"[SimNLLDerivativesHelper][init] --- LIST OF addition ---"< Date: Fri, 2 Jul 2021 18:23:54 +0200 Subject: [PATCH 09/16] work in progress --- interface/SimNLLDerivativesHelper.h | 2 + src/CascadeMinimizer.cc | 2 +- src/SimNLLDerivativesHelper.cc | 97 ++++++++++++++++++++++------- 3 files changed, 78 insertions(+), 23 deletions(-) diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h index ee0814a1231..8742ba2e608 100644 --- a/interface/SimNLLDerivativesHelper.h +++ b/interface/SimNLLDerivativesHelper.h @@ -57,6 +57,8 @@ class SimNLLDerivativesHelper //std::map& derivatives(){return derivatives;} private: cacheutils::CachingSimNLL* nll_{nullptr}; //not owned, keep pointer + + std::set getServersVars(RooAbsArg *node); }; #endif diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index 6c887fb1f8d..b100ae99fc8 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -79,7 +79,7 @@ void CascadeMinimizer::remakeMinimizer() { } SimNLLDerivativesHelper helper(simnll); helper.init(); - minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_,helper.derivatives_)); + minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_, (runtimedef::get("MINIMIZER_SemiAnalytic_NOANALYTIC"))?derivatives_:helper.derivatives_)); //minimizerSemiAnalytic_->setVerbose(kTRUE); // DEBUG } else{ diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index 03bf61d7e7b..a0651f270aa 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -141,13 +141,20 @@ void SimNLLDerivativesHelper::init(){ if ( dynamic_cast( &pp->components()[ip]) != nullptr) pn = dynamic_cast( &pp->components()[ip]); else { if (verbose) { std::cout <<"components:";} - RooArgList l = RooArgList(*pp->components()[ip].getVariables()); - for (int idx=0; idxGetName()) != logNormal.end()) { - logNormal.erase(v->GetName()); - if (verbose) std::cout<<"|"<GetName(); - } + //RooArgList l = RooArgList(*pp->components()[ip].getVariables()); + //for (int idx=0; idxGetName()) != logNormal.end()) { + // logNormal.erase(v->GetName()); + // if (verbose) std::cout<<"|"<GetName(); + // } + //} + std::set servers= getServersVars(&pp->components()[ip]); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|"<getVariables()); - for (int idx=0; idxGetName()) != logNormal.end()) { - logNormal.erase(v->GetName()); - if (verbose) std::cout<<"|"<GetName(); - } + //RooArgList l = RooArgList(*coeff->getVariables()); + //for (int idx=0; idxGetName()) != logNormal.end()) { + // logNormal.erase(v->GetName()); + // if (verbose) std::cout<<"|"<GetName(); + // } + //} + std::set servers= getServersVars(coeff); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|"<getVariables()); - for (int idx=0; idxGetName()) != logNormal.end()) { - logNormal.erase(v->GetName()); - if (verbose) std::cout<<"|"<GetName(); - } + + //l = RooArgList(*pdf->getVariables()); + //for (int idx=0; idxGetName()) != logNormal.end()) { + // logNormal.erase(v->GetName()); + // if (verbose) std::cout<<"|"<GetName(); + // } + //} + servers.clear(); + /*std::set*/ servers= getServersVars(pdf); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|"<> Removing all asymmThetaList "< SimNLLDerivativesHelper::getServersVars(RooAbsArg *node){ + std::set R; + auto iter = node->serverIterator(); + while (true) + { + auto server = iter->Next(); + if (server==nullptr)break; + + if (dynamic_cast (server) != nullptr) + { + R.insert(std::string(server->GetName())); + } + else{ + std::set r1 = getServersVars((RooAbsArg*)server); + R.merge(r1); + } + } + return R; +} From ffa5a4abf93b23f24f3bc8bbd1b77dd44bd6702e Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 6 Jul 2021 21:27:25 +0200 Subject: [PATCH 10/16] working for TH1 simple --- data/tutorials/shapes/simple-shapes-TH1.txt | 4 +- .../shapes/simple-shapes-parametric.txt | 2 +- interface/CascadeMinimizer.h | 2 + interface/RooMinimizerFcnSemiAnalytic.h | 6 +- interface/RooMinimizerSemiAnalytic.h | 2 +- interface/SimNLLDerivativesHelper.h | 44 +++-- interface/SimpleGaussianConstraint.h | 5 +- src/CachingNLL.cc | 2 +- src/CascadeMinimizer.cc | 7 +- src/RooMinimizerFcnSemiAnalytic.cxx | 58 ++++++- src/RooMinimizerSemiAnalytic.cc | 2 +- src/SimNLLDerivativesHelper.cc | 151 ++++++++++++------ 12 files changed, 201 insertions(+), 84 deletions(-) diff --git a/data/tutorials/shapes/simple-shapes-TH1.txt b/data/tutorials/shapes/simple-shapes-TH1.txt index 953eb92d017..c34f435af45 100644 --- a/data/tutorials/shapes/simple-shapes-TH1.txt +++ b/data/tutorials/shapes/simple-shapes-TH1.txt @@ -14,6 +14,6 @@ rate 10 100 -------------------------------- lumi lnN 1.10 1.0 bgnorm lnN 1.00 1.3 -alpha shapeN2 - 1 uncertainty on background shape and normalization -sigma shapeN2 0.5 - uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, +#alpha shapeN2 - 1 uncertainty on background shape and normalization +#sigma shapeN2 0.5 - uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, # so divide the unit gaussian by 2 before doing the interpolation diff --git a/data/tutorials/shapes/simple-shapes-parametric.txt b/data/tutorials/shapes/simple-shapes-parametric.txt index dc01d5e3320..6091d95cf0b 100644 --- a/data/tutorials/shapes/simple-shapes-parametric.txt +++ b/data/tutorials/shapes/simple-shapes-parametric.txt @@ -15,5 +15,5 @@ process sig bkg process 0 1 rate 1 1 -------------------------------- -lumi lnN 1.1 1.0 +lumi lnN 1.1 - vogian_sigma param 1.0 0.1 diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index e6575c7fd8e..3d3c61f3cce 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -13,6 +13,7 @@ class RooRealVar; #include #include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" #include class CascadeMinimizer { @@ -50,6 +51,7 @@ class CascadeMinimizer { bool isSemiAnalyticMinimizer{false}; std::auto_ptr minimizerSemiAnalytic_; std::map derivatives_; // not used + std::unique_ptr helper_; // Mode mode_; static int strategy_; diff --git a/interface/RooMinimizerFcnSemiAnalytic.h b/interface/RooMinimizerFcnSemiAnalytic.h index 16377ea3783..c73fdb9e000 100644 --- a/interface/RooMinimizerFcnSemiAnalytic.h +++ b/interface/RooMinimizerFcnSemiAnalytic.h @@ -37,7 +37,7 @@ class RooMinimizerFcnSemiAnalytic : public: - RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic *context, const std::map& knownDerivatives, bool verbose = false); + RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic *context, std::map* knownDerivatives, bool verbose = false); RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other); virtual ~RooMinimizerFcnSemiAnalytic(); @@ -90,7 +90,7 @@ class RooMinimizerFcnSemiAnalytic : RooAbsReal *_funct; // for name -> Derivative - std::map _knownDerivatives; // does not own them. Otherwise needs to design ad hoc copy constructors. + std::map *_knownDerivatives; // does not own them. Otherwise needs to design ad hoc copy constructors. RooMinimizerSemiAnalytic *_context; @@ -111,7 +111,7 @@ class RooMinimizerFcnSemiAnalytic : RooArgList* _initFloatParamList; RooArgList* _initConstParamList; - int _useNumDerivatives{1};// 0 = ROOT (not-impl), 1 re-implementation of gsl + int _useNumDerivatives{1};// 0 = ROOT (not-impl), 1 re-implementation of gsl 5 point, 2 re-implementation of gsl 3 point }; diff --git a/interface/RooMinimizerSemiAnalytic.h b/interface/RooMinimizerSemiAnalytic.h index 063d86447f4..d27df39c5dd 100644 --- a/interface/RooMinimizerSemiAnalytic.h +++ b/interface/RooMinimizerSemiAnalytic.h @@ -42,7 +42,7 @@ class RooPlot ; class RooMinimizerSemiAnalytic : public TObject { public: - RooMinimizerSemiAnalytic(RooAbsReal& function,const std::map& knownDerivatives) ; + RooMinimizerSemiAnalytic(RooAbsReal& function,std::map* knownDerivatives) ; virtual ~RooMinimizerSemiAnalytic() ; enum Strategy { Speed=0, Balance=1, Robustness=2 } ; diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h index 8742ba2e608..a7daedfdee2 100644 --- a/interface/SimNLLDerivativesHelper.h +++ b/interface/SimNLLDerivativesHelper.h @@ -14,8 +14,10 @@ #include #include #include +#include #include "RooAddition.h" #include "RooDataSet.h" +#include "RooRealVar.h" #include "HiggsAnalysis/CombinedLimit/interface/CachingNLL.h" #include "HiggsAnalysis/CombinedLimit/interface/ProcessNormalization.h" @@ -26,39 +28,59 @@ class DerivativeLogNormal: public RooAbsReal { const RooDataSet * data_{nullptr}; // not owned. - const cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL - std::string kappaname_{""}; + // RooRealProxy? + cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL + RooRealProxy pdfproxy_; + std::string thetaname_{""}; + + RooRealProxy theta_;// useful for numerical derivatives -> terms public: - DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& kappaname,int&found); + DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& thetaname,int&found); ~DerivativeLogNormal(){}; virtual Bool_t isDerived() const { return kTRUE; } virtual DerivativeLogNormal *clone(const char *name = 0) const ; const RooDataSet *data() const {return data_;} - const cacheutils::CachingAddNLL *pdf() const { return pdf_; } + cacheutils::CachingAddNLL *pdf() { return pdf_; } - virtual Double_t evaluate() const ; // cacheutils::ReminderSum + virtual Double_t evaluate() const override ; // cacheutils::ReminderSum // name-> index association per process std::vector kappa_pos_; bool verbose{true}; + + //bool maskConstraints{false}; + + // for debug purposes. Compute numerically the derivative corresponding to evaluate + //Double_t numericalDerivative() const; }; class SimNLLDerivativesHelper { + private: + //RooRealProxy? + cacheutils::CachingSimNLL* nll_{nullptr}; //not owned, keep pointer + RooRealVar unmask_;//("mask_derivative_constraint_","",1.); + + std::set getServersVars(const RooAbsArg *node); public: - SimNLLDerivativesHelper( cacheutils::CachingSimNLL * nll ){nll_=nll;} - ~SimNLLDerivativesHelper(){}; + SimNLLDerivativesHelper( cacheutils::CachingSimNLL * nll ): unmask_("mask_derivative_constraint_","",1.) {nll_=nll;} + ~SimNLLDerivativesHelper(){ + for (auto x : derivatives_) x.second->Delete(); + derivatives_.clear(); + }; void init(); bool verbose{true}; // debug std::map derivatives_; - //std::map& derivatives(){return derivatives;} - private: - cacheutils::CachingSimNLL* nll_{nullptr}; //not owned, keep pointer - std::set getServersVars(RooAbsArg *node); + void setMaskConstraint(int val=1){ + if (val ==0 or val ==1){ + unmask_.setVal( double(1-val) ); + } + else throw std::invalid_argument("[SimNLLDerivativesHelper]::setMask called with value != 0,1"); + } }; #endif diff --git a/interface/SimpleGaussianConstraint.h b/interface/SimpleGaussianConstraint.h index 364dd271c5a..a8bf9985b25 100644 --- a/interface/SimpleGaussianConstraint.h +++ b/interface/SimpleGaussianConstraint.h @@ -5,7 +5,7 @@ class SimpleGaussianConstraint : public RooGaussian { public: - SimpleGaussianConstraint() {} ; + SimpleGaussianConstraint() {}; SimpleGaussianConstraint(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma): RooGaussian(name,title,_x,_mean,_sigma) { init(); } @@ -26,7 +26,8 @@ class SimpleGaussianConstraint : public RooGaussian { //Double_t sig = sigma ; //return -0.5*arg*arg/(sig*sig); _value = scale_*arg*arg; - _valueDirty = false; + //_valueDirty = false; + _valueDirty = true; // something wrong with cache? } return _value; } diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index c8194eca339..0440a9017b4 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -1077,12 +1077,12 @@ cacheutils::CachingSimNLL::evaluate() const double logpdfval = (*it)->getLogValFast(); //std::cout << "pdf " << (*it)->GetName() << " = " << logpdfval << std::endl; ret2 += (logpdfval + *itz); + //std::cout << "[CachingSimNLL]: Simple Gaussian Constraint pdf " << (*it)->GetName() << " = " << logpdfval << " pars="<<(*it)->getMean().getVal()<<" "<<(*it)->getX().getVal()<<"/"<<(*it)->getSigma().getVal() <<" ZP "<<*itz<< std::endl; } /// ============= FAST POISSON CONSTRAINTS ========= itz = constrainZeroPointsFastPoisson_.begin(); for (std::vector::const_iterator it = constrainPdfsFastPoisson_.begin(), ed = constrainPdfsFastPoisson_.end(); it != ed; ++it, ++itz) { double logpdfval = (*it)->getLogValFast(); - //std::cout << "pdf " << (*it)->GetName() << " = " << logpdfval << std::endl; ret2 += (logpdfval + *itz); } } diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index b100ae99fc8..460dc7bd3dd 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -77,9 +77,10 @@ void CascadeMinimizer::remakeMinimizer() { std::cout<<"[ERROR]::[CascadeMinimizer] I need simnll to have SemiAnalytic derivatives done by SimNLLDerivativesHelper"<init(); + //helper_->setMaskConstraint(1); // Why the numerical derivatives correspond (wrongly) to this? TODO + minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_, (runtimedef::get("MINIMIZER_SemiAnalytic_NOANALYTIC"))? &derivatives_: &helper_->derivatives_)); //minimizerSemiAnalytic_->setVerbose(kTRUE); // DEBUG } else{ diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx index 12277dc2815..0c9bba770bb 100644 --- a/src/RooMinimizerFcnSemiAnalytic.cxx +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -36,11 +36,14 @@ #include "RooMinimizer.h" #include +#include "RooAddition.h" // DEBUG & Print +#include "RooFormulaVar.h" //DEBUG & Print +//#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" // DEBUG to have evaluate + using namespace std; -RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic* context, const std::map& knownDerivatives, bool verbose) : +RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic* context, std::map * knownDerivatives, bool verbose) : _funct(funct), - _knownDerivatives (knownDerivatives), _context(context), // Reset the *largest* negative log-likelihood value we have seen so far _maxFCN(-1e30), _numBadNLL(0), @@ -49,6 +52,7 @@ RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooM _verbose(verbose) { + _knownDerivatives = knownDerivatives; _evalCounter = 0 ; // Examine parameter list @@ -99,7 +103,6 @@ RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooM RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other) : ROOT::Math::IMultiGradFunction(other), _evalCounter(other._evalCounter), _funct(other._funct), - _knownDerivatives(other._knownDerivatives), _context(other._context), _maxFCN(other._maxFCN), _numBadNLL(other._numBadNLL), @@ -111,6 +114,7 @@ RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSe _floatParamVec(other._floatParamVec), _derivParamVec(other._derivParamVec) { + _knownDerivatives = other._knownDerivatives; _floatParamList = new RooArgList(*other._floatParamList) ; _constParamList = new RooArgList(*other._constParamList) ; _initFloatParamList = (RooArgList*) other._initFloatParamList->snapshot(kFALSE) ; @@ -520,8 +524,8 @@ void RooMinimizerFcnSemiAnalytic::updateFloatVec() _derivParamVec = std::vector(_floatParamList->getSize()) ; Int_t i(0) ; while((arg=iter.next())) { - auto derivative=_knownDerivatives.find(arg->GetName()); - if (derivative == _knownDerivatives.end()){ + auto derivative=_knownDerivatives->find(arg->GetName()); + if (derivative == _knownDerivatives->end()){ if(_verbose) std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "derivative for param "<GetName()<<" not found. I will use numerical derivatives."<getVal()<getVal()<getVal(); @@ -673,6 +682,7 @@ double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int i Compute the error using the difference between the 5-point and the 3-point rule (x-h,x,x+h). Again the central point is not used. */ + RooAbsReal::setHideOffset(kFALSE) ; // smallest number representable in doubles static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); @@ -698,6 +708,8 @@ double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int i SetPdfParamVal(icoord,x[icoord]+hhalf); double fph = _funct->getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + double r3 = 0.5 * (fp1 - fm1); double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3; @@ -717,8 +729,40 @@ double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int i //*result = r5 / h; //*abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */ //*abserr_round = fabs (e5 / h) + dy; /* Rounding error (cancellations) */ + // reset + SetPdfParamVal(icoord,x[icoord]); return r5 /h ; } + + if (_useNumDerivatives==2){ + /* Compute the derivative using the 5-point rule (x-h, x-h/2, x, + x+h/2, x+h). Note that the central point is not used. + + Compute the error using the difference between the 5-point and + the 3-point rule (x-h,x,x+h). Again the central point is not + used. */ + + // smallest number representable in doubles + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + static double fgEps = 0.001; + + double ax0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double h = (ax0>0) ? std::max( fgEps* ax0, 8.0*kPrecision*(ax0 + kPrecision) ) : kPrecision; + //static const double h=std::numeric_limits::epsilon()*8.; + + SetPdfParamVal(icoord,x[icoord]-h); + double fm1 = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+h); + double fp1 = _funct->getVal(); + + double r3 = 0.5 * (fp1 - fm1); + + // reset + SetPdfParamVal(icoord,x[icoord]); + return r3 /h ; + } throw std::runtime_error(Form("Numerical derivatives method un-implemented: %d",_useNumDerivatives)); } diff --git a/src/RooMinimizerSemiAnalytic.cc b/src/RooMinimizerSemiAnalytic.cc index 9b14684e24f..13e08fe3ba7 100644 --- a/src/RooMinimizerSemiAnalytic.cc +++ b/src/RooMinimizerSemiAnalytic.cc @@ -98,7 +98,7 @@ void RooMinimizerSemiAnalytic::cleanup() /// for HESSE and MINOS error analysis is taken from the defaultErrorLevel() /// value of the input function. -RooMinimizerSemiAnalytic::RooMinimizerSemiAnalytic(RooAbsReal& function,const std::map& knownDerivatives) +RooMinimizerSemiAnalytic::RooMinimizerSemiAnalytic(RooAbsReal& function, std::map* knownDerivatives) { RooSentinel::activate() ; diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index a0651f270aa..751044c51af 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -1,11 +1,16 @@ #include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" +#include "RooDataHist.h" +//#define DERIVATIVE_LOGNORMAL_DEBUG 1 -DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, const cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& kappaname,int&found) : +DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& thetaname,int&found) : RooAbsReal(name, title), data_(data), pdf_(pdf), - kappaname_(kappaname) + pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf), + thetaname_(thetaname) { + setDirtyInhibit(true); // TODO: understand cache and proxies + found=0; // keep track if at least one depends on the given kappa for (unsigned int i=0;ipdfs_.size();++i) // process loop { @@ -26,49 +31,76 @@ DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, co for (unsigned int j=0 ; jlogKappa_.size();++j) { - if ( p->thetaList_.at(j)->GetName() == kappaname_ ) {val=int(j);} + if ( p->thetaList_.at(j)->GetName() == thetaname_ ) {val=int(j);} } if (val==-1) { - if(verbose) std::cout<<"[DerivativeLogNormal]: Unable to find kappa corresponding to "<GetName() <GetName() <GetName() << "because not ProcessNormalization" <GetName() << "because not ProcessNormalization" <GetName() <GetName() < sum of expectations in the bin - * lambdat_b -> sum of expecations times logK + * lambdat_b -> sum of ( expecations times logK) */ double sum=0.0; // caching: TODO - for (int ib=0;ibnumEntries();++ib) +#ifdef DERIVATIVE_LOGNORMAL_DEBUG + if(verbose) { + std::cout<<"[DerivativeLogNormal]: data"<Print("V"); + std::cout<<"[DerivativeLogNormal]: --"< (data_); + + // FIX for ADDNLL_ROOREALSUM_KEEPZEROS + std::vector reminder(pdf_->pdfs_.size(),0.0); // the sums must do 1, the ib==ndata does that with data=0 + int ndata=data_->numEntries(); + + for (int ib=0;ib<=ndata;++ib) + //for (unsigned int ib=0;ibweights_.size();++ib) { - data_->get(ib) ; - double db = data_->weight(); - double lambda=pdf_->getVal(); + auto x= (ibget(ib): data_->get(0) ; + double db = (ibweight() : 0.; + double bw = (dh != nullptr) ? dh->binVolume(): 1.0; + + //double db = pdf_->weights_[ib]; + //double bw = pdf_->binWidths_[ib]; + + //double lambda=pdf_->getVal()*bw; // sum coeff*pdf double lambdat=0.0; + double lambda=0.;// debug: recompute lambda with loop. not sure that the term before include everything I want + + //bool recursive=true; + //double running_prod=1.0; for (unsigned int i=0;ipdfs_.size();++i) // loop over processes { + // find kappa ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); // fix RooProduct @@ -79,12 +111,40 @@ Double_t DerivativeLogNormal::evaluate() const { if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); } } - ///--- + ///--- from CachingAddNLL + //static bool expEventsNoNorm = runtimedef::get("ADDNLL_ROOREALSUM_NONORM"); + //--- double logK = (kappa_pos_[i]>=0)? c -> logKappa_ [kappa_pos_[i]] : 0.0; - lambdat+= pdf_->coeffs_[i]->getVal() * pdf_->pdfs_.at(i).pdf()->getVal() * logK; + if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" pdf is a RooRealSumPdf"<pdfs_.at(i).pdf()->getVal(x)*bw; + + double expectedEvents= (ibcoeffs_[i]->getVal(x) * pdf_->pdfs_.at(i).pdf()->getVal(x)*bw : + pdf_->coeffs_[i]->getVal(x) * (1.-reminder[i]); + + // sum coeff or sum coeff*integral + //double expectedEvents= dynamic_cast(pdf_->pdfs_.at(i).pdf())->expectedEvents(x); + //bool expEventsNoNorm=false; + //double expectedEvents = (pdf_->isRooRealSum_ && !expEventsNoNorm ? dynamic_cast(pdf_->pdfs_.at(i).pdf())->getNorm(data_->get()) : pdf_->coeffs_[i]->getVal()); + lambdat+= expectedEvents * logK; + lambda += expectedEvents; + + if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_< 1) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="<getVal()*bw <<")" <<" lambdat="< considering pdf "<GetName()<< ((isWeighted)?" isWeighted":" notWeighted") <pdfs_.size();++proc) { RooAbsReal * coeff = pdf -> coeffs_[proc] ; // what to do for RooProduct? + const RooAbsReal* pdfi= pdf->pdfs_[proc].pdf(); ProcessNormalization * pn = dynamic_cast (coeff); RooProduct *pp = dynamic_cast(coeff); @@ -141,14 +202,6 @@ void SimNLLDerivativesHelper::init(){ if ( dynamic_cast( &pp->components()[ip]) != nullptr) pn = dynamic_cast( &pp->components()[ip]); else { if (verbose) { std::cout <<"components:";} - //RooArgList l = RooArgList(*pp->components()[ip].getVariables()); - //for (int idx=0; idxGetName()) != logNormal.end()) { - // logNormal.erase(v->GetName()); - // if (verbose) std::cout<<"|"<GetName(); - // } - //} std::set servers= getServersVars(&pp->components()[ip]); for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { @@ -179,18 +232,12 @@ void SimNLLDerivativesHelper::init(){ if (verbose) { std::cout<<"[SimNLLDerivativesHelper][init]"<< " -- COEFF -- "<Print("V"); + std::cout<<"[SimNLLDerivativesHelper][init]"<< " -- PDF -- "<Print("V"); std::cout<<"-----------------"<getVariables()); - //for (int idx=0; idxGetName()) != logNormal.end()) { - // logNormal.erase(v->GetName()); - // if (verbose) std::cout<<"|"<GetName(); - // } - //} std::set servers= getServersVars(coeff); for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { @@ -201,16 +248,8 @@ void SimNLLDerivativesHelper::init(){ //-- if (verbose) { std::cout<getVariables()); - //for (int idx=0; idxGetName()) != logNormal.end()) { - // logNormal.erase(v->GetName()); - // if (verbose) std::cout<<"|"<GetName(); - // } - //} servers.clear(); - /*std::set*/ servers= getServersVars(pdf); + /*std::set*/ servers= getServersVars(pdfi); for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { logNormal.erase(v); @@ -222,9 +261,10 @@ void SimNLLDerivativesHelper::init(){ continue; // if pn==nullptr or not weighted } // not a process normalization coefficient + // TODO: remove all shapes + std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all asymmThetaList "<asymmThetaList_) for(int idx=0; idx< pn->asymmThetaList_.getSize() ;++idx){ RooAbsArg*v=pn->asymmThetaList_.at(idx); if (verbose) { @@ -235,7 +275,6 @@ void SimNLLDerivativesHelper::init(){ if (logNormal.find(v->GetName()) != logNormal.end()) logNormal.erase(v->GetName()); } // remove all otherFactorList - //for (auto v : pn->otherFactorList_) std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all otherfactorsList "<otherFactorList_.getSize() ;++idx){ RooAbsArg*v=pn->otherFactorList_.at(idx); @@ -272,13 +311,18 @@ void SimNLLDerivativesHelper::init(){ cacheutils::CachingAddNLL * pdf = *it; if (pdf==nullptr) continue ; // ?!? needed? const RooAbsData* data= pdf->data(); - //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string kappaname) + //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string thetaname) int found; - DerivativeLogNormal *der= new DerivativeLogNormal("","",pdf,dynamic_cast(data),name,found); - if(found) list.addOwned(*der); // or addOwned or add? + DerivativeLogNormal *der= new DerivativeLogNormal((std::string("dln_")+name+Form("_%d",idx)).c_str(),"",pdf,dynamic_cast(data),name,found); + if(found) list.add(*der); // or addOwned or add? else der->Delete(); } - // add derivative of constraint + + if (list.getSize() == 0){ // FIX for shapes TODO + std::cout<<"[SimNLLDerivativesHelper][init][ERROR]"<< "Something went wrong: Unable to match "<< name <<" with kappas"<constrainPdfsFast_ std::cout<<"[SimNLLDerivativesHelper][init]"<< "Adding constraint term for "<< name <constrainPdfsFast_.size();++i) @@ -287,8 +331,8 @@ void SimNLLDerivativesHelper::init(){ if (nll_->constrainPdfsFast_[i]->getX().GetName() != name) continue; if(verbose) std::cout<<"[SimNLLDerivativesHelper][init]"<< "Loop. found constraint "<< nll_->constrainPdfsFast_[i]->getX().GetName() <<" == "<constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma()));// (x-mean)/sigma^2 - list.addOwned(*constr) ; // or addOwned or add? + RooFormulaVar *constr = new RooFormulaVar( (std::string("constr_")+name).c_str(),"@3*(@0-@1)/(@2*@2)",RooArgList(nll_->constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma(),unmask_));// (x-mean)/sigma^2 + list.add(*constr) ; // or addOwned or add? } if (verbose) {std::cout <<"[SimNLLDerivativesHelper][init] --- LIST OF addition ---"<first <<" : "<second->getVal(); // DEBUG + } } //def getServers(node): @@ -312,7 +359,7 @@ void SimNLLDerivativesHelper::init(){ // servers.append(server) // return servers -std::set SimNLLDerivativesHelper::getServersVars(RooAbsArg *node){ +std::set SimNLLDerivativesHelper::getServersVars(const RooAbsArg *node){ std::set R; auto iter = node->serverIterator(); while (true) From f82e2643f91152a6738ce8b89fe7c176409ba0ee Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Tue, 6 Jul 2021 21:29:41 +0200 Subject: [PATCH 11/16] restore data/tutorials --- data/tutorials/shapes/simple-shapes-TH1.txt | 4 ++-- data/tutorials/shapes/simple-shapes-parametric.txt | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/data/tutorials/shapes/simple-shapes-TH1.txt b/data/tutorials/shapes/simple-shapes-TH1.txt index c34f435af45..953eb92d017 100644 --- a/data/tutorials/shapes/simple-shapes-TH1.txt +++ b/data/tutorials/shapes/simple-shapes-TH1.txt @@ -14,6 +14,6 @@ rate 10 100 -------------------------------- lumi lnN 1.10 1.0 bgnorm lnN 1.00 1.3 -#alpha shapeN2 - 1 uncertainty on background shape and normalization -#sigma shapeN2 0.5 - uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, +alpha shapeN2 - 1 uncertainty on background shape and normalization +sigma shapeN2 0.5 - uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, # so divide the unit gaussian by 2 before doing the interpolation diff --git a/data/tutorials/shapes/simple-shapes-parametric.txt b/data/tutorials/shapes/simple-shapes-parametric.txt index 6091d95cf0b..dc01d5e3320 100644 --- a/data/tutorials/shapes/simple-shapes-parametric.txt +++ b/data/tutorials/shapes/simple-shapes-parametric.txt @@ -15,5 +15,5 @@ process sig bkg process 0 1 rate 1 1 -------------------------------- -lumi lnN 1.1 - +lumi lnN 1.1 1.0 vogian_sigma param 1.0 0.1 From aa8dce68234dd4b3a9b6980eecdc946b1b41e511 Mon Sep 17 00:00:00 2001 From: Andrea Carlo Marini Date: Wed, 7 Jul 2021 11:19:20 +0200 Subject: [PATCH 12/16] fix for shapes uncertainties. still working on parametric cards. --- src/SimNLLDerivativesHelper.cc | 53 +++++++++++++++++++++++++--------- 1 file changed, 40 insertions(+), 13 deletions(-) diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index 751044c51af..9bd53c2fe7f 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -76,7 +76,8 @@ Double_t DerivativeLogNormal::evaluate() const { } #endif - const RooDataHist* dh = dynamic_cast (data_); + //const RooDataHist* dh = dynamic_cast (data_); + //if (dh==nullptr) std::cout<<"[DerivativeLogNormal]:"<<" NO ROODATAHIST"< reminder(pdf_->pdfs_.size(),0.0); // the sums must do 1, the ib==ndata does that with data=0 @@ -87,7 +88,23 @@ Double_t DerivativeLogNormal::evaluate() const { { auto x= (ibget(ib): data_->get(0) ; double db = (ibweight() : 0.; - double bw = (dh != nullptr) ? dh->binVolume(): 1.0; + //double bw = (dh != nullptr) ? dh->binVolume(): 1.0; + double bw = 1; //(pdf_->binWidths_.size() > 1) ? pdf_->binWidths_[ib]: (pdf_->binWidths_.size()==1) ?pdf_->binWidths_[0]: 1.0; + // const RooArgSet *obs = data_->get(); + // RooRealVar *xvar = dynamic_cast(obs->first()); + if (x->getSize()==1) { + RooRealVar *xvar=dynamic_cast(x->first()); + const RooAbsBinning &bins = xvar->getBinning(); + bw=bins.binWidth(0); // only costant supported. Need to figure out for zero bins + //if (verbose){ + // std::cout <<"BinWidth: "; + // for(int ii =0 ;iinumBins();++ii){ + // double bc2 = bins.binCenter(ii); + // std::cout<< bins.binWidth(ii); + // } + // std::cout <weights_[ib]; //double bw = pdf_->binWidths_[ib]; @@ -111,13 +128,12 @@ Double_t DerivativeLogNormal::evaluate() const { if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); } } - ///--- from CachingAddNLL - //static bool expEventsNoNorm = runtimedef::get("ADDNLL_ROOREALSUM_NONORM"); + //double integral = (pdf_->isRooRealSum_ && pdf_->basicIntegrals_ < 2) ? pdf_->integrals_[i]->getVal(x): 1.0; //--- double logK = (kappa_pos_[i]>=0)? c -> logKappa_ [kappa_pos_[i]] : 0.0; if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" pdf is a RooRealSumPdf"<pdfs_.at(i).pdf()->getVal(x)*bw; + if (ibpdfs_.at(i).pdf()->getVal(x)*bw; double expectedEvents= (ibcoeffs_[i]->getVal(x) * pdf_->pdfs_.at(i).pdf()->getVal(x)*bw : pdf_->coeffs_[i]->getVal(x) * (1.-reminder[i]); @@ -131,9 +147,9 @@ Double_t DerivativeLogNormal::evaluate() const { if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_< 1) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="< servers= getServersVars(pdfi); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|"<Print("V"); @@ -236,6 +264,7 @@ void SimNLLDerivativesHelper::init(){ pdfi->Print("V"); std::cout<<"-----------------"< servers= getServersVars(coeff); @@ -249,7 +278,7 @@ void SimNLLDerivativesHelper::init(){ if (verbose) { std::cout<*/ servers= getServersVars(pdfi); + /*std::set*/ servers= getServersVars(pdfi); // do I still need this block? for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { logNormal.erase(v); @@ -261,8 +290,6 @@ void SimNLLDerivativesHelper::init(){ continue; // if pn==nullptr or not weighted } // not a process normalization coefficient - // TODO: remove all shapes - std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all asymmThetaList "<asymmThetaList_.getSize() ;++idx){ @@ -318,7 +345,7 @@ void SimNLLDerivativesHelper::init(){ else der->Delete(); } - if (list.getSize() == 0){ // FIX for shapes TODO + if (list.getSize() == 0){ // old fix for shapes. Leave it, probably useful for debug std::cout<<"[SimNLLDerivativesHelper][init][ERROR]"<< "Something went wrong: Unable to match "<< name <<" with kappas"< Date: Fri, 16 Jul 2021 13:59:36 +0200 Subject: [PATCH 13/16] starting to implement derivatives for rate param (and simple pois) --- interface/CachingNLL.h | 2 + interface/ProcessNormalization.h | 2 + interface/SimNLLDerivativesHelper.h | 41 +++++++ src/RooMinimizerFcnSemiAnalytic.cxx | 9 +- src/SimNLLDerivativesHelper.cc | 165 ++++++++++++++++++++++++++-- 5 files changed, 207 insertions(+), 12 deletions(-) diff --git a/interface/CachingNLL.h b/interface/CachingNLL.h index 84e2806a0ab..43e89459253 100644 --- a/interface/CachingNLL.h +++ b/interface/CachingNLL.h @@ -23,6 +23,7 @@ class RooMultiPdf; class SimNLLDerivativesHelper; class DerivativeLogNormal; +class DerivativeRateParam; // Part zero: ArgSet checker namespace cacheutils { @@ -116,6 +117,7 @@ CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; class CachingAddNLL : public RooAbsReal { friend SimNLLDerivativesHelper; // probably not needed w/ data friend DerivativeLogNormal; + friend DerivativeRateParam; public: CachingAddNLL(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data, bool includeZeroWeights = false) ; CachingAddNLL(const CachingAddNLL &other, const char *name = 0) ; diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 61bf610f216..ca4cd7ab48f 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -16,6 +16,7 @@ END_HTML // class SimNLLDerivativesHelper; class DerivativeLogNormal; +class DerivativeRateParam; class ProcessNormalization : public RooAbsReal { public: @@ -38,6 +39,7 @@ class ProcessNormalization : public RooAbsReal { private: friend SimNLLDerivativesHelper; friend DerivativeLogNormal; + friend DerivativeRateParam; // ---- PERSISTENT ---- double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h index a7daedfdee2..e9939a9a121 100644 --- a/interface/SimNLLDerivativesHelper.h +++ b/interface/SimNLLDerivativesHelper.h @@ -25,6 +25,31 @@ * */ +// TODO: Make a common DerivativeAbstract class. Marke DerivativeLogNormal Inherits from this class +class DerivativeAbstract: public RooAbsReal +{ + protected: // Full access to derived classes + const RooDataSet * data_{nullptr}; // not owned. + // RooRealProxy? + cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL + RooRealProxy pdfproxy_; + public: + + DerivativeAbstract(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data): + RooAbsReal(name, title), + data_(data), + pdf_(pdf), + pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf) {}; + ~DerivativeAbstract(){}; + virtual Bool_t isDerived() const { return kTRUE; } + const RooDataSet *data() const {return data_;} + cacheutils::CachingAddNLL *pdf() { return pdf_; } + bool verbose{true}; + + virtual Double_t evaluate() const =0 ; // cacheutils::ReminderSum + +}; + class DerivativeLogNormal: public RooAbsReal { const RooDataSet * data_{nullptr}; // not owned. @@ -56,6 +81,22 @@ class DerivativeLogNormal: public RooAbsReal //Double_t numericalDerivative() const; }; +// Implementation of the derivative for rate params and rate POIs +class DerivativeRateParam : public DerivativeAbstract +{ + std::string ratename_{""}; + RooRealProxy rate_;// useful for numerical derivatives -> terms + public: + DerivativeRateParam(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& thetaname,int&found); + ~DerivativeRateParam(){}; + virtual DerivativeRateParam *clone(const char *name = 0) const ; + + virtual Double_t evaluate() const override ; // cacheutils::ReminderSum + // name-> index association per process in otherFactorList + std::vector rate_pos_; + +}; + class SimNLLDerivativesHelper { diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx index 0c9bba770bb..bc6a33c1370 100644 --- a/src/RooMinimizerFcnSemiAnalytic.cxx +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -39,6 +39,8 @@ #include "RooAddition.h" // DEBUG & Print #include "RooFormulaVar.h" //DEBUG & Print //#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" // DEBUG to have evaluate +#define DEBUG_DIRTY + using namespace std; @@ -609,6 +611,7 @@ double RooMinimizerFcnSemiAnalytic::DoEval(const double *x) const double RooMinimizerFcnSemiAnalytic::DoDerivative(const double * x, unsigned int icoord) const{ //std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "DoDerivative:"<getVal(); + + Double_t ret = _derivParamVec[icoord]->getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + return ret; } diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index 9bd53c2fe7f..3db3f6db8b8 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -1,6 +1,6 @@ #include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" #include "RooDataHist.h" -//#define DERIVATIVE_LOGNORMAL_DEBUG 1 +//#define DERIVATIVE_RATEPARAM_DEBUG 1 DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& thetaname,int&found) : RooAbsReal(name, title), @@ -9,7 +9,7 @@ DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, ca pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf), thetaname_(thetaname) { - setDirtyInhibit(true); // TODO: understand cache and proxies + //setDirtyInhibit(true); // TODO: understand cache and proxies found=0; // keep track if at least one depends on the given kappa for (unsigned int i=0;ipdfs_.size();++i) // process loop @@ -57,7 +57,7 @@ DerivativeLogNormal* DerivativeLogNormal::clone(const char *name) const { } Double_t DerivativeLogNormal::evaluate() const { -#ifdef DERIVATIVE_LOGNORMAL_DEBUG +#ifdef DERIVATIVE_RATEPARAM_DEBUG if(verbose) std::cout<<"[DerivativeLogNormal]: evaluate"<Print("V"); @@ -133,8 +133,10 @@ Double_t DerivativeLogNormal::evaluate() const { double logK = (kappa_pos_[i]>=0)? c -> logKappa_ [kappa_pos_[i]] : 0.0; if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" pdf is a RooRealSumPdf"<pdfs_.at(i).pdf()->getVal(x)*bw; + //pdf_->pdfs_.at(i).pdf()->setValueDirty(); + // TH1 + if (ibpdfs_.at(i).pdf()->getVal(x)*bw; double expectedEvents= (ibcoeffs_[i]->getVal(x) * pdf_->pdfs_.at(i).pdf()->getVal(x)*bw : pdf_->coeffs_[i]->getVal(x) * (1.-reminder[i]); @@ -147,18 +149,19 @@ Double_t DerivativeLogNormal::evaluate() const { if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_< 1*/) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<getVal()*bw <<")" <<" lambdat="<Print("V"); @@ -405,3 +408,145 @@ std::set SimNLLDerivativesHelper::getServersVars(const RooAbsArg * } return R; } + + + +/// --- +DerivativeRateParam::DerivativeRateParam(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& ratename,int&found) : + DerivativeAbstract(name,title,pdf,data), + ratename_(ratename) +{ + //setDirtyInhibit(true); // TODO: understand cache and proxies + + found=0; // keep track if at least one depends on the given kappa + for (unsigned int i=0;ipdfs_.size();++i) // process loop + { + if(verbose) std::cout<<"[DerivativeRateParam]: loop "<< i <<"/"<pdfs_.size()<<" considering pdf "<GetName() < (pdf_->coeffs_[i]); + if (verbose) std::cout<<"[DerivativeRateParam]: coeff is of class "<coeffs_[i]->ClassName()<(pdf_->coeffs_[i]); + if (p == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) p = dynamic_cast( &pp->components()[ip]); + } + } + // + if (p != nullptr) { // unable to add the sum for this process + + for (int j=0 ; jotherFactorList_.getSize();++j) + { + if ( p->otherFactorList_.at(j)->GetName() == ratename_ ) {val=int(j);} + } + if (val==-1) { + if(verbose) std::cout<<"[DerivativeRateParam]: Unable to find kappa corresponding to "<GetName() <GetName() << "because not ProcessNormalization" <GetName() < sum of expectations in the bin + * lambdat_b -> sum of ( expectations with delta(rate param)) + */ + + double sum=0.0; + // caching: TODO +#ifdef DERIVATIVE_RATEPARAM_DEBUG + if(verbose) { + std::cout<<"[DerivativeRateParam]: data"<Print("V"); + std::cout<<"[DerivativeRateParam]: --"< (data_); + //if (dh==nullptr) std::cout<<"[DerivativeRateParam]:"<<" NO ROODATAHIST"< reminder(pdf_->pdfs_.size(),0.0); // the sums must do 1, the ib==ndata does that with data=0 + int ndata=data_->numEntries(); + + for (int ib=0;ib<=ndata;++ib) + //for (unsigned int ib=0;ibweights_.size();++ib) + { + auto x= (ibget(ib): data_->get(0) ; + double db = (ibweight() : 0.; + double bw = 1; + if (x->getSize()==1) { + RooRealVar *xvar=dynamic_cast(x->first()); + const RooAbsBinning &bins = xvar->getBinning(); + bw=bins.binWidth(0); // only costant supported. Need to figure out for zero bins + } + + //double lambda=pdf_->getVal()*bw; // sum coeff*pdf + double lambdat=0.0; + double lambda=0.;// debug: recompute lambda with loop. not sure that the term before include everything I want + + for (unsigned int i=0;ipdfs_.size();++i) // loop over processes + { + // find rate param + ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); + + // fix RooProduct + RooProduct *pp = dynamic_cast(pdf_->coeffs_[i]); + if (c == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); + } + } + + //--- + double diracDelta = (rate_pos_[i]>=0)?1.0 : 0.0; + if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeRateParam]:DEBUG "<<" pdf is a RooRealSumPdf"<pdfs_.at(i).pdf()->getVal(x)*bw; + double expectedEvents= (ibcoeffs_[i]->getVal(x) * pdf_->pdfs_.at(i).pdf()->getVal(x)*bw : + pdf_->coeffs_[i]->getVal(x) * (1.-reminder[i]); + + lambdat+= expectedEvents * diracDelta; + lambda += expectedEvents; + + if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeRateParam]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_<getVal()*bw <<")" <<" lambdat="< Date: Mon, 19 Jul 2021 19:19:28 +0200 Subject: [PATCH 14/16] implementing rate parameter (v1) --- src/SimNLLDerivativesHelper.cc | 158 +++++++++++++++++++++++++++------ 1 file changed, 133 insertions(+), 25 deletions(-) diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc index 3db3f6db8b8..80dc5e6a32b 100644 --- a/src/SimNLLDerivativesHelper.cc +++ b/src/SimNLLDerivativesHelper.cc @@ -1,6 +1,8 @@ #include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" #include "RooDataHist.h" + //#define DERIVATIVE_RATEPARAM_DEBUG 1 +//#define DERIVATIVE_LOGNORMAL_DEBUG 1 DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& thetaname,int&found) : RooAbsReal(name, title), @@ -9,7 +11,7 @@ DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, ca pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf), thetaname_(thetaname) { - //setDirtyInhibit(true); // TODO: understand cache and proxies + setDirtyInhibit(true); // TODO: understand cache and proxies found=0; // keep track if at least one depends on the given kappa for (unsigned int i=0;ipdfs_.size();++i) // process loop @@ -57,7 +59,7 @@ DerivativeLogNormal* DerivativeLogNormal::clone(const char *name) const { } Double_t DerivativeLogNormal::evaluate() const { -#ifdef DERIVATIVE_RATEPARAM_DEBUG +#ifdef DERIVATIVE_LOGNORMAL_DEBUG if(verbose) std::cout<<"[DerivativeLogNormal]: evaluate"<Print("V"); @@ -149,19 +151,19 @@ Double_t DerivativeLogNormal::evaluate() const { if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_< 1*/) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<getVal()*bw <<")" <<" lambdat="< (coeff); RooProduct *pp = dynamic_cast(coeff); + RooRealVar *rv = dynamic_cast(coeff); + //if (rv !=nullptr) std::cout<<"[SimNLLDerivativesHelper]::[FIXME]"<<"derivative for "GetName()<<" is known, but discarded in the loop :( "<components().getSize(); ++ip) @@ -225,8 +253,12 @@ void SimNLLDerivativesHelper::init(){ for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { logNormal.erase(v); - if (verbose) std::cout<<"|"<Print("V"); @@ -274,8 +310,12 @@ void SimNLLDerivativesHelper::init(){ for(auto v : servers) { if (logNormal.find(v) != logNormal.end()) { logNormal.erase(v); - if (verbose) std::cout<<"|"<> Removing all asymmThetaList "<GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + if (rateParams.find(v->GetName()) != rateParams.end()) {rateParams.erase(v->GetName());} // probably not necessary, since this are constrained } + // remove all otherFactorList std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all otherfactorsList "<otherFactorList_.getSize() ;++idx){ RooAbsArg*v=pn->otherFactorList_.at(idx); if (verbose) { std::cout<<"[SimNLLDerivativesHelper][INFO]:" - <<"Removing "<GetName()<<" from the list of known derivatives because appears in an otherFactors " + <<"Removing "<GetName()<<" from the list of known derivatives (LN) because appears in an otherFactors " <GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + if (dynamic_cast(v)==nullptr and rateParams.find(v->GetName()) != rateParams.end()) { + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<< " : "<ClassName()<<" from the list of known derivatives (R) because appears in an otherFactors " + <GetName()); + } } + // remove all symmThetaList + for(int idx=0; idx< pn->thetaList_.getSize() ;++idx){ + RooAbsArg*v=pn->thetaList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName() + // << ","<GetX().GetName() + // << ","<GetMean().GetName() + // << ","<GetSigma().GetName() + <<" from the list of known derivatives because appears in an thetaList " + <GetName()) != rateParams.end()) {rateParams.erase(v->GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetX().GetName()) != rateParams.end()) {rateParams.erase(v->GetX().GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetMean().GetName()) != rateParams.end()) {rateParams.erase(v->GetMean().GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetSigma().GetName()) != rateParams.end()) {rateParams.erase(v->GetSigma().GetName());} // probably not necessary, since this are constrained + } std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Continue loop "< Continue loop "<::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + // construct the derivative as sum + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string thetaname) + int found; + DerivativeRateParam *der= new DerivativeRateParam((std::string("dln_")+name+Form("_%d",idx)).c_str(),"",pdf,dynamic_cast(data),name,found); + if(found) list.add(*der); // or addOwned or add? + else der->Delete(); + } + + if (list.getSize() == 0){ // old fix for shapes. Leave it, probably useful for debug + std::cout<<"[SimNLLDerivativesHelper][init][ERROR]"<< "Something went wrong: Unable to match "<< name <<" with kappas"<constrainPdfsFast_ + if (verbose) {std::cout <<"[SimNLLDerivativesHelper][init] --- LIST OF addition ---"< (pdf_->coeffs_[i]); + if (rv !=nullptr){ + val=-2; // the whole coeff is a RooRealVar + found=1; + if(verbose) std::cout<<"[DerivativeRateParam]: All Coeff corresponding to "<GetName() <<" is a RooRealVar"<otherFactorList_.getSize();++j) { if ( p->otherFactorList_.at(j)->GetName() == ratename_ ) {val=int(j);} } if (val==-1) { - if(verbose) std::cout<<"[DerivativeRateParam]: Unable to find kappa corresponding to "<GetName() <GetName() <GetName() << "because not ProcessNormalization" <GetName() << "because not ProcessNormalization" <GetName() <=0)?1.0 : 0.0; + double diracDelta = (rate_pos_[i]>=0 or rate_pos_[i]==-2)?1.0 : 0.0; if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeRateParam]:DEBUG "<<" pdf is a RooRealSumPdf"<getVal()*bw <<")" <<" lambdat="<getVal()*bw <<")" <<" lambdat="<pdfs_.size();++i) // process loop @@ -520,6 +521,7 @@ DerivativeRateParam::DerivativeRateParam(const char *name, const char *title, ca ratename_(ratename) { //setDirtyInhibit(true); // TODO: understand cache and proxies + setOperMode(RooAbsArg::ADirty,false); // Recursive, bool found=0; // keep track if at least one depends on the given kappa for (unsigned int i=0;ipdfs_.size();++i) // process loop @@ -609,8 +611,9 @@ Double_t DerivativeRateParam::evaluate() const { for (unsigned int i=0;ipdfs_.size();++i) // loop over processes { - // find rate param + // find rate param. TODO make it more efficient using rate_pos ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); + RooRealVar *rv = dynamic_cast(pdf_->coeffs_[i]); // fix RooProduct RooProduct *pp = dynamic_cast(pdf_->coeffs_[i]); @@ -620,17 +623,33 @@ Double_t DerivativeRateParam::evaluate() const { if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); } } + RooRealVar *rate = dynamic_cast( (rate_pos_[i]==-2)? rv: (rate_pos_[i]>=0) ? &c->otherFactorList_[rate_pos_[i]] : nullptr ); //--- - double diracDelta = (rate_pos_[i]>=0 or rate_pos_[i]==-2)?1.0 : 0.0; + int diracDelta = (rate_pos_[i]>=0 or rate_pos_[i]==-2)?1 : 0; if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeRateParam]:DEBUG "<<" pdf is a RooRealSumPdf"<basicIntegrals_<getVal()*bw <<")" <<" lambdat="<getVal()*bw <<")" <<" lambdat="<basicIntegrals_< 1*/) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<