diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index bd225ac37528e..9b428b3822e81 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -113,6 +113,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore RooGenFunction.h RooGenProdProj.h RooGlobalFunc.h + RooGradMinimizerFcn.h RooGrid.h RooHistError.h RooHistFunc.h @@ -333,6 +334,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/RooGenFunction.cxx src/RooGenProdProj.cxx src/RooGlobalFunc.cxx + src/RooGradMinimizerFcn.cxx src/RooGrid.cxx src/RooHist.cxx src/RooHistError.cxx diff --git a/roofit/roofitcore/inc/RooGradMinimizerFcn.h b/roofit/roofitcore/inc/RooGradMinimizerFcn.h new file mode 100644 index 0000000000000..829fc228d45d5 --- /dev/null +++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h @@ -0,0 +1,86 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * + * VC, Vince Croft, DIANA / NYU, vincent.croft@cern.ch * + * * + * * + * 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_GRAD_MINIMIZER_FCN +#define ROO_GRAD_MINIMIZER_FCN + +#include "RooAbsMinimizerFcn.h" + +#include "Minuit2/MnStrategy.h" +#include "Minuit2/MnMatrix.h" // MnAlgebraicVector +#include "Math/IFunction.h" // IMultiGradFunction +#include "Fit/ParameterSettings.h" +#include "Minuit2/NumericalDerivator.h" + +#include + +class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbsMinimizerFcn { +public: + RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose = false); + RooGradMinimizerFcn(const RooGradMinimizerFcn &other); + ROOT::Math::IMultiGradFunction *Clone() const override; + + ROOT::Minuit2::MnStrategy getStrategy() const; + double getErrorDef() const; + void setStrategy(int istrat); + + Bool_t Synchronize(std::vector ¶meter_settings, Bool_t optConst, + Bool_t verbose = kFALSE) override; + + void synchronizeGradientParameterSettings(std::vector ¶meter_settings) const; + + inline bool returnsInMinuit2ParameterSpace() const override { return true; } + inline unsigned int NDim() const override { return getNDim(); } + inline void setStepTolerance(double step_tolerance) const { _gradf.SetStepTolerance(step_tolerance); } + inline void setGradTolerance(double grad_tolerance) const { _gradf.SetGradTolerance(grad_tolerance); } + inline void setNcycles(unsigned int ncycles) const { _gradf.SetNCycles(ncycles); } + inline void setErrorLevel(double error_level) const { _gradf.SetErrorLevel(error_level); } + inline std::string getFunctionName() const override { return _funct->GetName(); } + inline std::string getFunctionTitle() const override { return _funct->GetTitle(); } + inline void setOffsetting(Bool_t flag) override { _funct->enableOffsetting(flag); } + +private: + void runDerivator(unsigned int i_component) const; + + void resetHasBeenCalculatedFlags() const; + bool syncParameter(double x, std::size_t ix) const; + bool syncParameters(const double *x) const; + + inline void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override + { + _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); + } + +public: + enum class GradientCalculatorMode { ExactlyMinuit2, AlmostMinuit2 }; + +private: + // IMultiGradFunction overrides + double DoEval(const double *x) const override; + double DoDerivative(const double *x, unsigned int icoord) const override; + + // members + // mutable because ROOT::Math::IMultiGradFunction::DoDerivative is const +protected: + mutable std::vector _grad; + mutable std::vector _grad_params; + +private: + mutable ROOT::Minuit2::NumericalDerivator _gradf; + RooAbsReal *_funct; + mutable std::vector has_been_calculated; + mutable bool none_have_been_calculated = false; +}; +#endif diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 01c377fadeb06..835e71b0aa791 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -6,6 +6,7 @@ * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * * * * * * Redistribution and use in source and binary forms, * @@ -16,6 +17,8 @@ #ifndef ROO_MINIMIZER #define ROO_MINIMIZER +#include // shared_ptr, unique_ptr + #include "TObject.h" #include "TStopwatch.h" #include @@ -24,8 +27,16 @@ #include #include "TMatrixDSymfwd.h" -#include "Fit/Fitter.h" +#include "RooArgList.h" // cannot just use forward decl due to default argument in lastMinuitFit + #include "RooMinimizerFcn.h" +#include "RooGradMinimizerFcn.h" + +#include "RooSentinel.h" +#include "RooMsgService.h" + +#include "Fit/Fitter.h" +#include // logic_error class RooAbsReal ; class RooFitResult ; @@ -37,9 +48,11 @@ class RooPlot ; class RooMinimizer : public TObject { public: + enum class FcnMode { classic, gradient }; - RooMinimizer(RooAbsReal& function) ; - virtual ~RooMinimizer() ; + explicit RooMinimizer(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic); + static std::unique_ptr create(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic); + ~RooMinimizer() override; enum Strategy { Speed=0, Balance=1, Robustness=2 } ; enum PrintLevel { None=-1, Reduced=0, Normal=1, ExtraForProblem=2, Maximum=3 } ; @@ -91,7 +104,10 @@ class RooMinimizer : public TObject { ROOT::Fit::Fitter* fitter() ; const ROOT::Fit::Fitter* fitter() const ; - + + ROOT::Math::IMultiGenFunction* getFitterMultiGenFcn() const; + ROOT::Math::IMultiGenFunction* getMultiGenFcn() const; + protected: friend class RooAbsPdf ; @@ -100,28 +116,30 @@ class RooMinimizer : public TObject { void profileStart() ; void profileStop() ; - inline Int_t getNPar() const { return fitterFcn()->NDim() ; } + inline Int_t getNPar() const { return fitterFcn()->getNDim() ; } inline std::ofstream* logfile() { return fitterFcn()->GetLogFile(); } inline Double_t& maxFCN() { return fitterFcn()->GetMaxFCN() ; } - - const RooMinimizerFcn* fitterFcn() const { return ( fitter()->GetFCN() ? ((RooMinimizerFcn*) fitter()->GetFCN()) : _fcn ) ; } - RooMinimizerFcn* fitterFcn() { return ( fitter()->GetFCN() ? ((RooMinimizerFcn*) fitter()->GetFCN()) : _fcn ) ; } -private: + const RooAbsMinimizerFcn *fitterFcn() const; + RooAbsMinimizerFcn *fitterFcn(); - Int_t _printLevel ; - Int_t _status ; - Bool_t _profile ; + bool fitFcn() const; - Bool_t _verbose ; - TStopwatch _timer ; - TStopwatch _cumulTimer ; - Bool_t _profileStart ; +private: + Int_t _printLevel = 1; + Int_t _status = -99; + Bool_t _profile = kFALSE; + + Bool_t _verbose = kFALSE; + TStopwatch _timer; + TStopwatch _cumulTimer; + Bool_t _profileStart = kFALSE; - TMatrixDSym* _extV ; + TMatrixDSym *_extV = 0; - RooMinimizerFcn *_fcn; - std::string _minimizerType; + RooAbsMinimizerFcn *_fcn; + std::string _minimizerType = "Minuit"; + FcnMode _fcnMode; static ROOT::Fit::Fitter *_theFitter ; @@ -129,8 +147,7 @@ class RooMinimizer : public TObject { RooMinimizer(const RooMinimizer&) ; - ClassDef(RooMinimizer,0) // RooFit interface to ROOT::Fit::Fitter + ClassDefOverride(RooMinimizer,0) // RooFit interface to ROOT::Fit::Fitter } ; - #endif diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx new file mode 100644 index 0000000000000..55e5d66ab5010 --- /dev/null +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -0,0 +1,238 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * + * VC, Vince Croft, DIANA / NYU, vincent.croft@cern.ch * + * * + * * + * 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) * + *****************************************************************************/ + +////////////////////////////////////////////////////////////////////////////// +// +// RooGradMinimizerFcn is am interface class to the ROOT::Math function +// for minimization. See RooGradMinimizer.cxx for more information. +// + +#include "RooGradMinimizerFcn.h" + +#include "RooAbsArg.h" +#include "RooAbsPdf.h" +#include "RooArgSet.h" +#include "RooRealVar.h" +#include "RooMsgService.h" +#include "RooMinimizer.h" + +#include "Riostream.h" +#include "TIterator.h" +#include "TClass.h" +#include "Fit/Fitter.h" +#include "Math/Minimizer.h" + +#include // std::equal +#include + +RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose) + : RooAbsMinimizerFcn(RooArgList( * std::unique_ptr(funct->getParameters(RooArgSet{})) ), context, verbose), + _grad(getNDim()), _grad_params(getNDim()), _funct(funct), + has_been_calculated(getNDim()) +{ + // TODO: added "parameters" after rewrite in april 2020, check if correct + auto parameters = _context->fitter()->Config().ParamsSettings(); + synchronizeParameterSettings(parameters, kTRUE, verbose); + synchronizeGradientParameterSettings(parameters); + setStrategy(ROOT::Math::MinimizerOptions::DefaultStrategy()); + setErrorLevel(ROOT::Math::MinimizerOptions::DefaultErrorDef()); +} + +RooGradMinimizerFcn::RooGradMinimizerFcn(const RooGradMinimizerFcn &other) + : RooAbsMinimizerFcn(other), _grad(other._grad), _grad_params(other._grad_params), _gradf(other._gradf), _funct(other._funct), + has_been_calculated(other.has_been_calculated), none_have_been_calculated(other.none_have_been_calculated) +{ +} + +ROOT::Math::IMultiGradFunction *RooGradMinimizerFcn::Clone() const +{ + return new RooGradMinimizerFcn(*this); +} + +void RooGradMinimizerFcn::synchronizeGradientParameterSettings( + std::vector ¶meter_settings) const +{ + _gradf.SetInitialGradient(_context->getMultiGenFcn(), parameter_settings, _grad); +} + +//////////////////////////////////////////////////////////////////////////////// + +double RooGradMinimizerFcn::DoEval(const double *x) const +{ + Bool_t parameters_changed = kFALSE; + + // Set the parameter values for this iteration + for (unsigned index = 0; index < NDim(); index++) { + // also check whether the function was already evaluated for this set of parameters + parameters_changed |= SetPdfParamVal(index, x[index]); + } + + // Calculate the function for these parameters + RooAbsReal::setHideOffset(kFALSE); + double fvalue = _funct->getVal(); + RooAbsReal::setHideOffset(kTRUE); + + if (!parameters_changed) { + return fvalue; + } + + if (!std::isfinite(fvalue) || RooAbsReal::numEvalErrors() > 0 || fvalue > 1e30) { + + if (_printEvalErrors >= 0) { + + if (_doEvalErrorWall) { + oocoutW(static_cast(nullptr), Eval) + << "RooGradMinimizerFcn: Minimized function has error status." << std::endl + << "Returning maximum FCN so far (" << _maxFCN + << ") to force MIGRAD to back out of this region. Error log follows" << std::endl; + } else { + oocoutW(static_cast(nullptr), Eval) + << "RooGradMinimizerFcn: Minimized function has error status but is ignored" << std::endl; + } + + TIterator *iter = _floatParamList->createIterator(); + RooRealVar *var; + Bool_t first(kTRUE); + ooccoutW(static_cast(nullptr), Eval) << "Parameter values: "; + while ((var = (RooRealVar *)iter->Next())) { + if (first) { + first = kFALSE; + } else + ooccoutW(static_cast(nullptr), Eval) << ", "; + ooccoutW(static_cast(nullptr), Eval) << var->GetName() << "=" << var->getVal(); + } + delete iter; + ooccoutW(static_cast(nullptr), Eval) << std::endl; + + RooAbsReal::printEvalErrors(ooccoutW(static_cast(nullptr), Eval), _printEvalErrors); + ooccoutW(static_cast(nullptr), Eval) << std::endl; + } + + if (_doEvalErrorWall) { + fvalue = _maxFCN + 1; + } + + RooAbsReal::clearEvalErrorLog(); + _numBadNLL++; + } else if (fvalue > _maxFCN) { + _maxFCN = fvalue; + } + + // Optional logging + if (_verbose) { + std::cout << "\nprevFCN" << (_funct->isOffsetting() ? "-offset" : "") << " = " << std::setprecision(10) << fvalue + << std::setprecision(4) << " "; + std::cout.flush(); + } + + _evalCounter++; + //#ifndef NDEBUG + // std::cout << "RooGradMinimizerFcn " << this << " evaluations (in DoEval): " << _evalCounter << + // std::endl; + //#endif + return fvalue; +} + +void RooGradMinimizerFcn::resetHasBeenCalculatedFlags() const +{ + for (auto it = has_been_calculated.begin(); it != has_been_calculated.end(); ++it) { + *it = false; + } + none_have_been_calculated = true; +} + +bool RooGradMinimizerFcn::syncParameter(double x, std::size_t ix) const +{ + bool parameter_has_changed = (_grad_params[ix] != x); + + if (parameter_has_changed) { + _grad_params[ix] = x; + // Set the parameter values for this iteration + // TODO: this is already done in DoEval as well; find efficient way to do only once + SetPdfParamVal(ix, x); + + if (!none_have_been_calculated) { + resetHasBeenCalculatedFlags(); + } + } + + return parameter_has_changed; +} + +bool RooGradMinimizerFcn::syncParameters(const double *x) const +{ + bool has_been_synced = false; + + for (std::size_t ix = 0; ix < NDim(); ++ix) { + bool parameter_has_changed = (_grad_params[ix] != x[ix]); + + if (parameter_has_changed) { + _grad_params[ix] = x[ix]; + // Set the parameter values for this iteration + // TODO: this is already done in DoEval as well; find efficient way to do only once + SetPdfParamVal(ix, x[ix]); + } + + has_been_synced |= parameter_has_changed; + } + + if (has_been_synced) { + resetHasBeenCalculatedFlags(); + } + + return has_been_synced; +} + +void RooGradMinimizerFcn::runDerivator(unsigned int i_component) const +{ + // check whether the derivative was already calculated for this set of parameters + if (!has_been_calculated[i_component]) { + // Calculate the derivative etc for these parameters + _grad[i_component] = + _gradf.PartialDerivative(_context->getMultiGenFcn(), _grad_params.data(), + _context->fitter()->Config().ParamsSettings(), i_component, _grad[i_component]); + has_been_calculated[i_component] = true; + none_have_been_calculated = false; + } +} + +double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int i_component) const +{ + syncParameters(x); + runDerivator(i_component); + return _grad[i_component].derivative; +} + +//////////////////////////////////////////////////////////////////////////////// + +void RooGradMinimizerFcn::setStrategy(int istrat) +{ + assert(istrat >= 0); + ROOT::Minuit2::MnStrategy strategy(static_cast(istrat)); + + setStepTolerance(strategy.GradientStepTolerance()); + setGradTolerance(strategy.GradientTolerance()); + setNcycles(strategy.GradientNCycles()); +} + +Bool_t +RooGradMinimizerFcn::Synchronize(std::vector ¶meters, Bool_t optConst, Bool_t verbose) +{ + Bool_t returnee = synchronizeParameterSettings(parameters, optConst, verbose); + synchronizeGradientParameterSettings(parameters); + setStrategy(_context->fitter()->Config().MinimizerOptions().Strategy()); + setErrorLevel(_context->fitter()->Config().MinimizerOptions().ErrorDef()); + return returnee; +} diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index a579d85bbdd8b..7bcdab4c708c7 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -6,6 +6,7 @@ * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * * * * Redistribution and use in source and binary forms, * * with or without modification, are permitted according to the terms * @@ -20,7 +21,8 @@ RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between the minimizer functionality and the native RooFit interface. -By default the Minimizer is MINUIT. +By default the Minimizer is MINUIT for classic FcnMode and MINUIT2 +for gradient FcnMode. RooMinimizer can minimize any RooAbsReal function with respect to its parameters. Usual choices for minimization are RooNLLVar and RooChi2Var @@ -36,6 +38,7 @@ Various methods are available to control verbosity, profiling, automatic PDF optimization. **/ + #include "RooFit.h" #include "TClass.h" @@ -60,9 +63,8 @@ automatic PDF optimization. #include "RooMsgService.h" #include "RooPlot.h" - -#include "RooMinimizer.h" #include "RooFitResult.h" +#include "RooMinimizer.h" #include "Math/Minimizer.h" @@ -105,45 +107,57 @@ void RooMinimizer::cleanup() /// for HESSE and MINOS error analysis is taken from the defaultErrorLevel() /// value of the input function. -RooMinimizer::RooMinimizer(RooAbsReal& function) +RooMinimizer::RooMinimizer(RooAbsReal &function, FcnMode fcnMode) : _fcnMode(fcnMode) { - RooSentinel::activate() ; - - // Store function reference - _extV = 0 ; - _verbose = kFALSE ; - _profile = kFALSE ; - _profileStart = kFALSE ; - _printLevel = 1 ; - _minimizerType = "Minuit"; // default minimizer - - if (_theFitter) delete _theFitter ; - _theFitter = new ROOT::Fit::Fitter; - _fcn = new RooMinimizerFcn(&function,this,_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(function.defaultErrorLevel()) ; - - // Declare our parameters to MINUIT - _fcn->Synchronize(_theFitter->Config().ParamsSettings(), - _fcn->getOptConst(),_verbose) ; - - // Now set default verbosity - if (RooMsgService::instance().silentMode()) { - setPrintLevel(-1) ; - } else { - setPrintLevel(1) ; - } + RooSentinel::activate(); + + if (_theFitter) { + delete _theFitter; + } + _theFitter = new ROOT::Fit::Fitter; + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + setEps(1.0); // default tolerance + + switch (_fcnMode) { + case FcnMode::classic: { + _fcn = new RooMinimizerFcn(&function, this, _verbose); + break; + } + case FcnMode::gradient: { + _fcn = new RooGradMinimizerFcn(&function, this, _verbose); + setMinimizerType("Minuit2"); + break; + } + default: { + throw std::logic_error("In RooMinimizer constructor: fcnMode has an unsupported value!"); + } + } + + // default max number of calls + _theFitter->Config().MinimizerOptions().SetMaxIterations(500 * _fcn->getNDim()); + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500 * _fcn->getNDim()); + + // Shut up for now + setPrintLevel(-1); + + // Use +0.5 for 1-sigma errors + setErrorLevel(function.defaultErrorLevel()); + + // Declare our parameters to MINUIT + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), _fcn->getOptConst(), _verbose); + + // Now set default verbosity + if (RooMsgService::instance().silentMode()) { + setPrintLevel(-1); + } else { + setPrintLevel(1); + } } +// static function +std::unique_ptr RooMinimizer::create(RooAbsReal &function, FcnMode fcnMode) { + return std::make_unique(function, fcnMode); +} //////////////////////////////////////////////////////////////////////////////// @@ -234,12 +248,17 @@ void RooMinimizer::setOffsetting(Bool_t flag) } +//////////////////////////////////////////////////////////////////////////////// +/// Choose the minimizer algorithm. +// forward declaration (necessary for avoiding circular dependency problems) +class RooGradMinimizerFcn; -//////////////////////////////////////////////////////////////////////////////// -/// Choose the minimiser algorithm. void RooMinimizer::setMinimizerType(const char* type) { + if (_fcnMode != FcnMode::classic && strcmp(type, "Minuit2") != 0) { + throw std::invalid_argument("In RooMinimizer::setMinimizerType: only Minuit2 is supported when not using classic function mode!"); + } _minimizerType = type; } @@ -298,6 +317,25 @@ RooFitResult* RooMinimizer::fit(const char* options) } +bool RooMinimizer::fitFcn() const { + bool ret; + + switch (_fcnMode) { + case FcnMode::classic: { + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); + break; + } + case FcnMode::gradient: { + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); + break; + } + default: { + throw std::logic_error("In RooMinimizer::fitFcn: _fcnMode has an unsupported value!"); + } + } + + return ret; +} //////////////////////////////////////////////////////////////////////////////// @@ -317,7 +355,7 @@ Int_t RooMinimizer::minimize(const char* type, const char* alg) RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; RooAbsReal::clearEvalErrorLog() ; - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = fitFcn(); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -346,7 +384,7 @@ Int_t RooMinimizer::migrad() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = fitFcn(); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -509,7 +547,7 @@ Int_t RooMinimizer::seek() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = fitFcn(); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -538,7 +576,7 @@ Int_t RooMinimizer::simplex() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = fitFcn(); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -567,7 +605,7 @@ Int_t RooMinimizer::improve() RooAbsReal::clearEvalErrorLog() ; _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved"); - bool ret = _theFitter->FitFCN(*_fcn); + bool ret = fitFcn(); _status = ((ret) ? _theFitter->Result().Status() : -1); RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; @@ -816,7 +854,56 @@ void RooMinimizer::profileStop() } +ROOT::Math::IMultiGenFunction* RooMinimizer::getFitterMultiGenFcn() const +{ + return fitter()->GetFCN(); +} + + +ROOT::Math::IMultiGenFunction* RooMinimizer::getMultiGenFcn() const +{ + if (getFitterMultiGenFcn()) { + return getFitterMultiGenFcn(); + } else { + switch (_fcnMode) { + case FcnMode::classic: { + return static_cast(dynamic_cast(_fcn)); + } + case FcnMode::gradient: { + return static_cast(dynamic_cast(_fcn)); + } + default: { + throw std::logic_error("In RooMinimizer::getMultiGenFcn: _fcnMode has an unsupported value!"); + } + } + } +} + +const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const +{ + if (getFitterMultiGenFcn()) { + switch (_fcnMode) { + case FcnMode::classic: { + return static_cast(dynamic_cast(getFitterMultiGenFcn())); + } + case FcnMode::gradient: { + return static_cast(dynamic_cast(getFitterMultiGenFcn())); + } + default: { + throw std::logic_error("In RooMinimizer::fitterFcn: _fcnMode has an unsupported value!"); + } + } + } else { + return _fcn; + } +} + +RooAbsMinimizerFcn *RooMinimizer::fitterFcn() +{ + // to avoid code duplication, we just reuse the const function and cast constness away + return const_cast( static_cast(*this).fitterFcn() ); +} //////////////////////////////////////////////////////////////////////////////// diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 9ec4dd851a82a..9a1a54c81d4d1 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -5,6 +5,7 @@ # For the list of contributors see $ROOTSYS/README/CREDITS. # @author Danilo Piparo CERN, 2018 +# @author Patrick Bos, NL eScience Center, 2018 ROOT_ADD_GTEST(simple simple.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooCacheManager testRooCacheManager.cxx LIBRARIES RooFitCore) @@ -39,4 +40,5 @@ endif() ROOT_ADD_GTEST(testRooProductPdf testRooProductPdf.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore RooFit) +#ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore) diff --git a/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx new file mode 100644 index 0000000000000..3d1a8181612c8 --- /dev/null +++ b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx @@ -0,0 +1,647 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * + * * + * 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) * + *****************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include "TFile.h" // for loading the workspace file +#include // remove redundant workspace files + +#include + +#include "gtest/gtest.h" +#include "test_lib.h" + +#include +#include // RooFit::ERROR + +TEST(GradMinimizer, Gaussian1D) +{ + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + for (int i = 0; i < 10; ++i) { + RooWorkspace w = RooWorkspace(); + + std::unique_ptr nll; + std::unique_ptr values; + RooAbsPdf *pdf; + RooDataSet *data; + std::tie(nll, pdf, data, values) = generate_1D_gaussian_pdf_nll(w, 10000); + // when c++17 support arrives, change to this: + // auto [nll, values] = generate_1D_gaussian_pdf_nll(w, 10000); + RooRealVar *mu = w.var("mu"); + + RooArgSet *savedValues = dynamic_cast(values->snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // -------- + + RooMinimizer m0(*nll); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(-1); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mu0 = mu->getVal(); + double muerr0 = mu->getError(); + + *values = *savedValues; + + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(-1); + + m1->migrad(); + + RooFitResult *m1result = m1->lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mu1 = mu->getVal(); + double muerr1 = mu->getError(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(mu0, mu1); + EXPECT_EQ(muerr0, muerr1); + EXPECT_EQ(edm0, edm1); + } +} + +TEST(GradMinimizerDebugging, DISABLED_Gaussian1DNominal) +{ + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w = RooWorkspace(); + + std::unique_ptr nll; + std::unique_ptr _; + RooAbsPdf *pdf; + RooDataSet *data; + std::tie(nll, pdf, data, _) = generate_1D_gaussian_pdf_nll(w, 10000); + // when c++17 support arrives, change to this: + // auto [nll, _] = generate_1D_gaussian_pdf_nll(w, 10000); + + RooMinimizer m0(*nll); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(100); + m0.setVerbose(kTRUE); + + m0.migrad(); +} + +TEST(GradMinimizerDebugging, DISABLED_Gaussian1DGradMinimizer) +{ + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w = RooWorkspace(); + + std::unique_ptr nll; + std::unique_ptr _; + RooAbsPdf *pdf; + RooDataSet *data; + std::tie(nll, pdf, data, _) = generate_1D_gaussian_pdf_nll(w, 10000); + // when c++17 support arrives, change to this: + // auto [nll, _] = generate_1D_gaussian_pdf_nll(w, 10000); + + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(100); + m1->setVerbose(kTRUE); + + m1->migrad(); +} + +TEST(GradMinimizer, GaussianND) +{ + // test RooMinimizer class with simple N-dimensional pdf + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + unsigned N = 5; + unsigned N_events = 1000; + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w("w", kFALSE); + + std::unique_ptr nll; + std::unique_ptr all_values; + RooAbsPdf *pdf; + RooDataSet *data; + std::tie(nll, pdf, data, all_values) = generate_ND_gaussian_pdf_nll(w, N, N_events); + // when c++17 support arrives, change to this: + // auto [nll, all_values] = generate_ND_gaussian_pdf_nll(w, N, N_events); + + // save initial values for the start of all minimizations + RooArgSet *savedValues = dynamic_cast(all_values->snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // -------- + + RooMinimizer m0(*(nll.get())); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(-1); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mean0[N]; + double std0[N]; + for (unsigned ix = 0; ix < N; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + mean0[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + { + std::ostringstream os; + os << "s" << ix; + std0[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + } + + // -------- + + *all_values = *savedValues; + + // -------- + + std::unique_ptr m1 = RooMinimizer::create(*(nll.get()), RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(-1); + + m1->migrad(); + + RooFitResult *m1result = m1->lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mean1[N]; + double std1[N]; + for (unsigned ix = 0; ix < N; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + mean1[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + { + std::ostringstream os; + os << "s" << ix; + std1[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + } + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(edm0, edm1); + + for (unsigned ix = 0; ix < N; ++ix) { + EXPECT_EQ(mean0[ix], mean1[ix]); + EXPECT_EQ(std0[ix], std1[ix]); + } +} + +TEST(GradMinimizerReverse, GaussianND) +{ + // test RooMinimizer class with simple N-dimensional pdf + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + unsigned N = 5; + unsigned N_events = 1000; + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w("w", kFALSE); + + std::unique_ptr nll; + std::unique_ptr all_values; + RooAbsPdf *pdf; + RooDataSet *data; + std::tie(nll, pdf, data, all_values) = generate_ND_gaussian_pdf_nll(w, N, N_events); + // when c++17 support arrives, change to this: + // auto [nll, all_values] = generate_ND_gaussian_pdf_nll(w, N, N_events); + + // save initial values for the start of all minimizations + RooArgSet *savedValues = dynamic_cast(all_values->snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // -------- + + std::unique_ptr m0 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + + m0->setMinimizerType("Minuit2"); + + m0->setStrategy(0); + m0->setPrintLevel(-1); + + m0->migrad(); + + RooFitResult *m0result = m0->lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mean0[N]; + double std0[N]; + for (unsigned ix = 0; ix < N; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + mean0[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + { + std::ostringstream os; + os << "s" << ix; + std0[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + } + + // -------- + + *all_values = *savedValues; + + // -------- + + RooMinimizer m1(*nll); + m1.setMinimizerType("Minuit2"); + + m1.setStrategy(0); + m1.setPrintLevel(-1); + + m1.migrad(); + + RooFitResult *m1result = m1.lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mean1[N]; + double std1[N]; + for (unsigned ix = 0; ix < N; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + mean1[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + { + std::ostringstream os; + os << "s" << ix; + std1[ix] = dynamic_cast(w.arg(os.str().c_str()))->getVal(); + } + } + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(edm0, edm1); + + for (unsigned ix = 0; ix < N; ++ix) { + EXPECT_EQ(mean0[ix], mean1[ix]); + EXPECT_EQ(std0[ix], std1[ix]); + } +} + +TEST(GradMinimizer, BranchingPDF) +{ + // test RooMinimizer class with an N-dimensional pdf that forms a tree of + // pdfs, where one subpdf is the parameter of a higher level pdf + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + int N_events = 1000; + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w("w", kFALSE); + + // 3rd level + w.factory("Gamma::ga0_0_1(k0_0_1[3,2,10],u[1,20],1,0)"); // leaf pdf + // Gamma(mu,N+1,1,0) ~ Pois(N,mu), so this is a "continuous Poissonian" + + // 2nd level that will be linked to from 3rd level + w.factory("Gamma::ga1_0(k1_0[4,2,10],z[1,20],1,0)"); // leaf pdf + + // rest of 3rd level + w.factory("Gaussian::g0_0_0(v[-10,10],m0_0_0[0.6,-10,10],ga1_0)"); // two branch pdf, one up a level to different 1st + // level branch + + // rest of 2nd level + w.factory("Gaussian::g0_0(g0_0_0,m0_0[6,-10,10],ga0_0_1)"); // branch pdf + + // 1st level + w.factory("Gaussian::g0(x[-10,10],g0_0,s0[3,0.1,10])"); // branch pdf + w.factory("Gaussian::g1(y[-10,10],m1[-2,-10,10],ga1_0)"); // branch pdf + RooArgSet level1_pdfs; + level1_pdfs.add(*w.arg("g0")); + level1_pdfs.add(*w.arg("g1")); + + // event counts for 1st level pdfs + RooRealVar N_g0("N_g0", "#events g0", N_events / 10, 0., 10 * N_events); + RooRealVar N_g1("N_g1", "#events g1", N_events / 10, 0., 10 * N_events); + w.import(N_g0); + w.import(N_g1); + // gather in count_set + RooArgSet level1_counts; + level1_counts.add(N_g0); + level1_counts.add(N_g1); + + // finally, sum the top level pdfs + RooAddPdf sum("sum", "gaussian tree", level1_pdfs, level1_counts); + + // gather observables + RooArgSet obs_set; + for (auto obs : {"x", "y", "z", "u", "v"}) { + obs_set.add(*w.arg(obs)); + } + + // --- Generate a toyMC sample from composite PDF --- + RooDataSet *data = sum.generate(obs_set, N_events); + + auto nll = sum.createNLL(*data); + + // gather all values of parameters, observables, pdfs and nll here for easy + // saving and restoring + RooArgSet some_values = RooArgSet(obs_set, w.allPdfs(), "some_values"); + RooArgSet most_values = RooArgSet(some_values, level1_counts, "most_values"); + most_values.add(*nll); + most_values.add(sum); + + RooArgSet *param_set = nll->getParameters(obs_set); + + RooArgSet all_values = RooArgSet(most_values, *param_set, "all_values"); + + // set parameter values randomly so that they actually need to do some fitting + auto it = all_values.fwdIterator(); + while (RooRealVar *val = dynamic_cast(it.next())) { + val->setVal(RooRandom::randomGenerator()->Uniform(val->getMin(), val->getMax())); + } + + // save initial values for the start of all minimizations + RooArgSet *savedValues = dynamic_cast(all_values.snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // write the workspace state to file, but remove it again if everything was + // successful (at the end of the test) + w.import(*data); + w.import(sum); + w.writeToFile("failed_testRooGradMinimizer_BranchingPDF_workspace.root"); + + // -------- + + RooMinimizer m0(*nll); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(-1); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + + double N_g0__0 = N_g0.getVal(); + double N_g1__0 = N_g1.getVal(); + double k0_0_1__0 = dynamic_cast(w.arg("k0_0_1"))->getVal(); + double k1_0__0 = dynamic_cast(w.arg("k1_0"))->getVal(); + double m0_0__0 = dynamic_cast(w.arg("m0_0"))->getVal(); + double m0_0_0__0 = dynamic_cast(w.arg("m0_0_0"))->getVal(); + double m1__0 = dynamic_cast(w.arg("m1"))->getVal(); + double s0__0 = dynamic_cast(w.arg("s0"))->getVal(); + + // -------- + + all_values = *savedValues; + + // -------- + + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(-1); + + m1->migrad(); + + RooFitResult *m1result = m1->lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(edm0, edm1); + + double N_g0__1 = N_g0.getVal(); + double N_g1__1 = N_g1.getVal(); + double k0_0_1__1 = dynamic_cast(w.arg("k0_0_1"))->getVal(); + double k1_0__1 = dynamic_cast(w.arg("k1_0"))->getVal(); + double m0_0__1 = dynamic_cast(w.arg("m0_0"))->getVal(); + double m0_0_0__1 = dynamic_cast(w.arg("m0_0_0"))->getVal(); + double m1__1 = dynamic_cast(w.arg("m1"))->getVal(); + double s0__1 = dynamic_cast(w.arg("s0"))->getVal(); + + EXPECT_EQ(N_g0__0, N_g0__1); + EXPECT_EQ(N_g1__0, N_g1__1); + EXPECT_EQ(k0_0_1__0, k0_0_1__1); + EXPECT_EQ(k1_0__0, k1_0__1); + EXPECT_EQ(m0_0__0, m0_0__1); + EXPECT_EQ(m0_0_0__0, m0_0_0__1); + EXPECT_EQ(m1__0, m1__1); + EXPECT_EQ(s0__0, s0__1); + + // N_g0 = 494.514 +/- 18.8621 (limited) + // N_g1 = 505.817 +/- 24.6705 (limited) + // k0_0_1 = 2.96883 +/- 0.00561152 (limited) + // k1_0 = 4.12068 +/- 0.0565994 (limited) + // m0_0 = 8.09563 +/- 1.30395 (limited) + // m0_0_0 = 0.411472 +/- 0.183239 (limited) + // m1 = -1.99988 +/- 0.00194089 (limited) + // s0 = 3.04623 +/- 0.0982477 (limited) + + if (!HasFailure()) { + if (remove("failed_testRooGradMinimizer_BranchingPDF_workspace.root") != 0) { + std::cout << "Failed to remove failed_testRooGradMinimizer_BranchingPDF_workspace.root workspace file, sorry. " + "There were no failures though, so manually remove at your leisure." + << std::endl; + } + } +} + +TEST(GradMinimizerDebugging, DISABLED_BranchingPDFLoadFromWorkspace) +{ + // test RooMinimizer class with an N-dimensional pdf that forms a tree of + // pdfs, where one subpdf is the parameter of a higher level pdf + + // This version of the BranchingPDF test loads the random parameters written + // to a workspace file by the original BranchingPDF test at some point. + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + TFile *f = new TFile("failed_testRooGradMinimizer_BranchingPDF_workspace.root"); + + // Retrieve workspace from file + RooWorkspace w = *static_cast(f->Get("w")); + + RooAddPdf sum = *static_cast(w.pdf("sum")); + RooDataSet *data = static_cast(w.data("")); + + auto nll = sum.createNLL(*data); + + RooArgSet all_values = w.allVars(); + // RooArgSet all_values; + // for (auto var_name : {"x", "y", "z", "u", "v", "ga0_0_1", "ga1_0", "g0_0_0", "g0_0", "g0", "g1", "N_g0", "N_g1", + // "nll_sum_sumData", "sum", "k0_0_1", "k1_0", "m0_0", "m0_0_0", "m1", "s0"}) { + // all_values.add(*w.arg(var_name)); + // } + + // save initial values for the start of all minimizations + RooArgSet *savedValues = dynamic_cast(all_values.snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // -------- + + all_values.Print("v"); + RooMinimizer m0(*nll); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(-1); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + + double N_g0__0 = dynamic_cast(w.arg("N_g0"))->getVal(); + double N_g1__0 = dynamic_cast(w.arg("N_g1"))->getVal(); + double k0_0_1__0 = dynamic_cast(w.arg("k0_0_1"))->getVal(); + double k1_0__0 = dynamic_cast(w.arg("k1_0"))->getVal(); + double m0_0__0 = dynamic_cast(w.arg("m0_0"))->getVal(); + double m0_0_0__0 = dynamic_cast(w.arg("m0_0_0"))->getVal(); + double m1__0 = dynamic_cast(w.arg("m1"))->getVal(); + double s0__0 = dynamic_cast(w.arg("s0"))->getVal(); + + all_values.Print("v"); + + // -------- + + all_values = *savedValues; + + // -------- + + all_values.Print("v"); + + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(-1); + + m1->migrad(); + + RooFitResult *m1result = m1->lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(edm0, edm1); + + double N_g0__1 = dynamic_cast(w.arg("N_g0"))->getVal(); + double N_g1__1 = dynamic_cast(w.arg("N_g1"))->getVal(); + double k0_0_1__1 = dynamic_cast(w.arg("k0_0_1"))->getVal(); + double k1_0__1 = dynamic_cast(w.arg("k1_0"))->getVal(); + double m0_0__1 = dynamic_cast(w.arg("m0_0"))->getVal(); + double m0_0_0__1 = dynamic_cast(w.arg("m0_0_0"))->getVal(); + double m1__1 = dynamic_cast(w.arg("m1"))->getVal(); + double s0__1 = dynamic_cast(w.arg("s0"))->getVal(); + + EXPECT_EQ(N_g0__0, N_g0__1); + EXPECT_EQ(N_g1__0, N_g1__1); + EXPECT_EQ(k0_0_1__0, k0_0_1__1); + EXPECT_EQ(k1_0__0, k1_0__1); + EXPECT_EQ(m0_0__0, m0_0__1); + EXPECT_EQ(m0_0_0__0, m0_0_0__1); + EXPECT_EQ(m1__0, m1__1); + EXPECT_EQ(s0__0, s0__1); + + all_values.Print("v"); + + // N_g0 = 494.514 +/- 18.8621 (limited) + // N_g1 = 505.817 +/- 24.6705 (limited) + // k0_0_1 = 2.96883 +/- 0.00561152 (limited) + // k1_0 = 4.12068 +/- 0.0565994 (limited) + // m0_0 = 8.09563 +/- 1.30395 (limited) + // m0_0_0 = 0.411472 +/- 0.183239 (limited) + // m1 = -1.99988 +/- 0.00194089 (limited) + // s0 = 3.04623 +/- 0.0982477 (limited) +} + +TEST(GradMinimizerDebugging, DISABLED_BranchingPDFLoadFromWorkspaceNominal) +{ + // only run the nominal minimizer of the BranchingPDF test and print results + + TFile *f = new TFile("failed_testRooGradMinimizer_BranchingPDF_workspace.root"); + RooWorkspace w = *static_cast(f->Get("w")); + RooAddPdf sum = *static_cast(w.pdf("sum")); + RooDataSet *data = static_cast(w.data("")); + auto nll = sum.createNLL(*data); + + RooMinimizer m0(*nll); + m0.setMinimizerType("Minuit2"); + m0.setStrategy(0); + m0.migrad(); +} + +TEST(GradMinimizerDebugging, DISABLED_BranchingPDFLoadFromWorkspaceGradMinimizer) +{ + // only run the GradMinimizer from the BranchingPDF test and print results + + TFile *f = new TFile("failed_testRooGradMinimizer_BranchingPDF_workspace.root"); + RooWorkspace w = *static_cast(f->Get("w")); + RooAddPdf sum = *static_cast(w.pdf("sum")); + RooDataSet *data = static_cast(w.data("")); + auto nll = sum.createNLL(*data); + + std::unique_ptr m0 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m0->setMinimizerType("Minuit2"); + m0->setStrategy(0); + m0->migrad(); +} diff --git a/roofit/roofitcore/test/test_lib.h b/roofit/roofitcore/test/test_lib.h new file mode 100644 index 0000000000000..dae7ef3b50dd7 --- /dev/null +++ b/roofit/roofitcore/test/test_lib.h @@ -0,0 +1,166 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * PB, Patrick Bos, NL eScience Center, p.bos@esciencecenter.nl * + * * + * 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 ROOT_TEST_LIB_H +#define ROOT_TEST_LIB_H + +#include "RooWorkspace.h" +#include "RooRandom.h" +#include "RooAddPdf.h" +#include "RooDataSet.h" +#include "RooRealVar.h" // for the dynamic cast to have a complete type + +#include +#include // make_unique + +RooAbsPdf * generate_1D_gaussian_pdf(RooWorkspace &w) +{ + w.factory("Gaussian::g(x[-5,5],mu[0,-3,3],sigma[1])"); + RooAbsPdf *pdf = w.pdf("g"); + return pdf; +} + +RooDataSet * generate_1D_dataset(RooWorkspace &w, RooAbsPdf *pdf, unsigned long N_events) +{ + RooDataSet *data = pdf->generate(RooArgSet(*w.var("x")), N_events); + return data; +} + + +std::tuple, RooAbsPdf *, RooDataSet *, std::unique_ptr> +generate_1D_gaussian_pdf_nll(RooWorkspace &w, unsigned long N_events) +{ + RooAbsPdf *pdf = generate_1D_gaussian_pdf(w); + + RooDataSet *data = generate_1D_dataset(w, pdf, N_events); + + RooRealVar *mu = w.var("mu"); + mu->setVal(-2.9); + + std::unique_ptr nll{pdf->createNLL(*data)}; + + // save initial values for the start of all minimizations + std::unique_ptr values = std::make_unique(*mu, *pdf, *nll, "values"); + + return std::make_tuple(std::move(nll), pdf, data, std::move(values)); +} + +// return two unique_ptrs, the first because nll is a pointer, +// the second because RooArgSet doesn't have a move ctor +std::tuple, RooAbsPdf *, RooDataSet *, std::unique_ptr> +generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_events) { + RooArgSet obs_set; + + // create gaussian parameters + double mean[n], sigma[n]; + for (unsigned ix = 0; ix < n; ++ix) { + mean[ix] = RooRandom::randomGenerator()->Gaus(0, 2); + sigma[ix] = 0.1 + abs(RooRandom::randomGenerator()->Gaus(0, 2)); + } + + // create gaussians and also the observables and parameters they depend on + for (unsigned ix = 0; ix < n; ++ix) { + std::ostringstream os; + os << "Gaussian::g" << ix + << "(x" << ix << "[-10,10]," + << "m" << ix << "[" << mean[ix] << ",-10,10]," + << "s" << ix << "[" << sigma[ix] << ",0.1,10])"; + w.factory(os.str().c_str()); + } + + // create uniform background signals on each observable + for (unsigned ix = 0; ix < n; ++ix) { + { + std::ostringstream os; + os << "Uniform::u" << ix << "(x" << ix << ")"; + w.factory(os.str().c_str()); + } + + // gather the observables in a list for data generation below + { + std::ostringstream os; + os << "x" << ix; + obs_set.add(*w.arg(os.str().c_str())); + } + } + + RooArgSet pdf_set = w.allPdfs(); + + // create event counts for all pdfs + RooArgSet count_set; + + // ... for the gaussians + for (unsigned ix = 0; ix < n; ++ix) { + std::stringstream os, os2; + os << "Nsig" << ix; + os2 << "#signal events comp " << ix; + RooRealVar a(os.str().c_str(), os2.str().c_str(), N_events/10, 0., 10*N_events); + w.import(a); + // gather in count_set + count_set.add(*w.arg(os.str().c_str())); + } + // ... and for the uniform background components + for (unsigned ix = 0; ix < n; ++ix) { + std::stringstream os, os2; + os << "Nbkg" << ix; + os2 << "#background events comp " << ix; + RooRealVar a(os.str().c_str(), os2.str().c_str(), N_events/10, 0., 10*N_events); + w.import(a); + // gather in count_set + count_set.add(*w.arg(os.str().c_str())); + } + + RooAddPdf* sum = new RooAddPdf("sum", "gaussians+uniforms", pdf_set, count_set); + w.import(*sum); // keep sum around after returning + + // --- Generate a toyMC sample from composite PDF --- + RooDataSet *data = sum->generate(obs_set, N_events); + + std::unique_ptr nll {sum->createNLL(*data)}; + + // set values randomly so that they actually need to do some fitting + for (unsigned ix = 0; ix < n; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + dynamic_cast(w.arg(os.str().c_str()))->setVal(RooRandom::randomGenerator()->Gaus(0, 2)); + } + { + std::ostringstream os; + os << "s" << ix; + dynamic_cast(w.arg(os.str().c_str()))->setVal(0.1 + abs(RooRandom::randomGenerator()->Gaus(0, 2))); + } + } + + // gather all values of parameters, pdfs and nll here for easy + // saving and restoring + std::unique_ptr all_values = std::make_unique(pdf_set, count_set, "all_values"); + all_values->add(*nll); + all_values->add(*sum); + for (unsigned ix = 0; ix < n; ++ix) { + { + std::ostringstream os; + os << "m" << ix; + all_values->add(*w.arg(os.str().c_str())); + } + { + std::ostringstream os; + os << "s" << ix; + all_values->add(*w.arg(os.str().c_str())); + } + } + + return std::make_tuple(std::move(nll), sum, data, std::move(all_values)); +} + + +#endif //ROOT_TEST_LIB_H