From f80ad0fad2520fdac68cb12624f804fb1cb07193 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Sat, 3 Jul 2021 21:13:47 +0200 Subject: [PATCH 1/7] add RooGradMinimizerFcn This commit adds RooGradMinimizerFcn, an alternative to RooMinimizerFcn that calculates gradient itself outside of Minuit during minimization. To be able to use RooGradMinimizerFcn, RooMinimizer had to be refactored to be able to use RooAbsMinimizerFcn instead of the concrete RooMinimizerFcn class. To accomodate this change, several other changes had to be made: - Functions getMultiGenFcn, fitterFcn and fitFcn were added which are used in several places as convenience functions to access the concrete RooAbsMinimizerFcn object with its proper type. The concrete type matters, for instance, when passing the class to the derivator, because a gradient enabled *MinimizerFcn needs to take different overloads than the classic function-value-only RooMinimizerFcn. A FcnMode enum class was added for this disambiguation as well. - Construction of a RooMinimizer can be done the same as always when using a classic RooMinimizerFcn. However, to activate RooGradMinimizerFcn, a new create function can be used with the proper type as template parameter. In a later PR, this pattern will be extended with another gradient-enabled RooAbsMinimizerFcn instantiation that will also include multi-processing implementations. Since we had to refactor construction, we took the opportunity to also put default parameters in the header declaration, instead of in the constructor initializer lists. Finally, this commit changes the default minimizer from Minuit to Minuit2. Other than these RooMinimizer changes and the addition of RooGradMinimizerFcn, a new test was added called testRooGradMinimizer. It also uses a test_lib.h header that will also be used later on in tests of other new likelihood minimization implementations. --- roofit/roofitcore/CMakeLists.txt | 2 + roofit/roofitcore/inc/RooGradMinimizerFcn.h | 92 ++ roofit/roofitcore/inc/RooMinimizer.h | 54 +- roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 294 ++++++ roofit/roofitcore/src/RooMinimizer.cxx | 152 ++-- roofit/roofitcore/test/CMakeLists.txt | 2 + roofit/roofitcore/test/RooGradMinimizer.cxx | 847 ++++++++++++++++++ roofit/roofitcore/test/test_lib.h | 185 ++++ 8 files changed, 1554 insertions(+), 74 deletions(-) create mode 100644 roofit/roofitcore/inc/RooGradMinimizerFcn.h create mode 100644 roofit/roofitcore/src/RooGradMinimizerFcn.cxx create mode 100644 roofit/roofitcore/test/RooGradMinimizer.cxx create mode 100644 roofit/roofitcore/test/test_lib.h 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..87605092de2fd --- /dev/null +++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h @@ -0,0 +1,92 @@ +/***************************************************************************** + * 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 + +#include "Minuit2/MnStrategy.h" +#include "Minuit2/MnMatrix.h" // MnAlgebraicVector +#include "Math/IFunction.h" // IMultiGradFunction +#include "Fit/ParameterSettings.h" +#include "Minuit2/NumericalDerivator.h" + +#include "RooAbsMinimizerFcn.h" + +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 get_strategy() const; + double get_error_def() const; + void set_strategy(int istrat); + + Bool_t Synchronize(std::vector ¶meter_settings, Bool_t optConst, + Bool_t verbose = kFALSE) override; + + void synchronize_gradient_parameter_settings(std::vector ¶meter_settings) const; + + bool returnsInMinuit2ParameterSpace() const override; + + unsigned int NDim() const override; + + void set_step_tolerance(double step_tolerance) const; + void set_grad_tolerance(double grad_tolerance) const; + void set_ncycles(unsigned int ncycles) const; + void set_error_level(double error_level) const; + + std::string getFunctionName() const override; + std::string getFunctionTitle() const override; + + void setOffsetting(Bool_t flag) override; + +private: + void run_derivator(unsigned int i_component) const; + + void reset_has_been_calculated_flags() const; + bool sync_parameter(double x, std::size_t ix) const; + bool sync_parameters(const double *x) const; + + void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override; + +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; + +public: + // for debugging, wraps ROOT::Math::IMultiGradFunction::Gradient, can be used for further inspection: + void Gradient(const double *x, double *grad) const override; +}; +#endif diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 01c377fadeb06..6eb930aa733a7 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,8 +48,11 @@ class RooPlot ; class RooMinimizer : public TObject { public: + enum class FcnMode { classic, gradient }; RooMinimizer(RooAbsReal& function) ; + template + static std::unique_ptr create(RooAbsReal &function); virtual ~RooMinimizer() ; enum Strategy { Speed=0, Balance=1, Robustness=2 } ; @@ -91,7 +105,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 +117,33 @@ 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 ) ; } + + const RooAbsMinimizerFcn *fitterFcn() const; + RooAbsMinimizerFcn *fitterFcn(); + + bool fitFcn() const; private: + template + RooMinimizer(RooAbsReal &function, MinimizerFcn* /* used only for template deduction */); - Int_t _printLevel ; - Int_t _status ; - Bool_t _profile ; + Int_t _printLevel = 1; + Int_t _status = -99; + Bool_t _profile = kFALSE; - Bool_t _verbose ; - TStopwatch _timer ; - TStopwatch _cumulTimer ; - Bool_t _profileStart ; + 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 = "Minuit2"; + FcnMode _fcnMode; static ROOT::Fit::Fitter *_theFitter ; diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx new file mode 100644 index 0000000000000..08c259d186810 --- /dev/null +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -0,0 +1,294 @@ +/***************************************************************************** + * 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 + +#include "Riostream.h" + +#include "TIterator.h" +#include "TClass.h" + +#include "RooAbsArg.h" +#include "RooAbsPdf.h" +#include "RooArgSet.h" +#include "RooRealVar.h" +#include "RooMsgService.h" + +#include "RooMinimizer.h" + +#include "Fit/Fitter.h" +#include "Math/Minimizer.h" + +#include // std::equal + +RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose) + : RooAbsMinimizerFcn(RooArgList(*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); + synchronize_gradient_parameter_settings(parameters); + set_strategy(ROOT::Math::MinimizerOptions::DefaultStrategy()); + set_error_level(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::synchronize_gradient_parameter_settings( + 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::reset_has_been_calculated_flags() const +{ + for (auto it = has_been_calculated.begin(); it != has_been_calculated.end(); ++it) { + *it = false; + } + none_have_been_calculated = true; +} + +bool RooGradMinimizerFcn::sync_parameter(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) { + reset_has_been_calculated_flags(); + } + } + + return parameter_has_changed; +} + +bool RooGradMinimizerFcn::sync_parameters(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) { + reset_has_been_calculated_flags(); + } + + return has_been_synced; +} + +void RooGradMinimizerFcn::run_derivator(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.partial_derivative(_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 +{ + sync_parameters(x); + run_derivator(i_component); + return _grad[i_component].derivative; +} + +bool RooGradMinimizerFcn::returnsInMinuit2ParameterSpace() const +{ + return true; +} + +unsigned int RooGradMinimizerFcn::NDim() const +{ + return getNDim(); +} + +//////////////////////////////////////////////////////////////////////////////// + +void RooGradMinimizerFcn::set_strategy(int istrat) +{ + assert(istrat >= 0); + ROOT::Minuit2::MnStrategy strategy(static_cast(istrat)); + + set_step_tolerance(strategy.GradientStepTolerance()); + set_grad_tolerance(strategy.GradientTolerance()); + set_ncycles(strategy.GradientNCycles()); +} + +Bool_t +RooGradMinimizerFcn::Synchronize(std::vector ¶meters, Bool_t optConst, Bool_t verbose) +{ + Bool_t returnee = synchronizeParameterSettings(parameters, optConst, verbose); + synchronize_gradient_parameter_settings(parameters); + set_strategy(_context->fitter()->Config().MinimizerOptions().Strategy()); + set_error_level(_context->fitter()->Config().MinimizerOptions().ErrorDef()); + return returnee; +} + +void RooGradMinimizerFcn::setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) +{ + _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); +} + +void RooGradMinimizerFcn::set_step_tolerance(double step_tolerance) const +{ + _gradf.set_step_tolerance(step_tolerance); +} +void RooGradMinimizerFcn::set_grad_tolerance(double grad_tolerance) const +{ + _gradf.set_grad_tolerance(grad_tolerance); +} +void RooGradMinimizerFcn::set_ncycles(unsigned int ncycles) const +{ + _gradf.set_ncycles(ncycles); +} +void RooGradMinimizerFcn::set_error_level(double error_level) const +{ + _gradf.set_error_level(error_level); +} + +std::string RooGradMinimizerFcn::getFunctionName() const +{ + return _funct->GetName(); +} + +std::string RooGradMinimizerFcn::getFunctionTitle() const +{ + return _funct->GetTitle(); +} + +void RooGradMinimizerFcn::setOffsetting(Bool_t flag) +{ + _funct->enableOffsetting(flag); +} + +void RooGradMinimizerFcn::Gradient(const double *x, double *grad) const +{ + ROOT::Math::IMultiGradFunction::Gradient(x, grad); +} diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index a579d85bbdd8b..df1bd6d8f1316 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,7 @@ 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 MINUIT2. RooMinimizer can minimize any RooAbsReal function with respect to its parameters. Usual choices for minimization are RooNLLVar and RooChi2Var @@ -36,6 +37,7 @@ Various methods are available to control verbosity, profiling, automatic PDF optimization. **/ + #include "RooFit.h" #include "TClass.h" @@ -60,9 +62,8 @@ automatic PDF optimization. #include "RooMsgService.h" #include "RooPlot.h" - -#include "RooMinimizer.h" #include "RooFitResult.h" +#include "RooMinimizer.h" #include "Math/Minimizer.h" @@ -94,58 +95,14 @@ void RooMinimizer::cleanup() //////////////////////////////////////////////////////////////////////////////// -/// 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. - -RooMinimizer::RooMinimizer(RooAbsReal& function) -{ - 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) ; +/// Constructor - // Now set default verbosity - if (RooMsgService::instance().silentMode()) { - setPrintLevel(-1) ; - } else { - setPrintLevel(1) ; - } +RooMinimizer::RooMinimizer(RooAbsReal &function) : + RooMinimizer::RooMinimizer(function, static_cast(nullptr)) +{ } - //////////////////////////////////////////////////////////////////////////////// /// Destructor @@ -234,12 +191,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 (dynamic_cast(_fcn) && strcmp(type, "Minuit2") != 0) { + throw std::invalid_argument("In RooMinimizer::setMinimizerType: only Minuit2 is supported when using RooGradMinimizerFcn!"); + } _minimizerType = type; } @@ -298,6 +260,27 @@ RooFitResult* RooMinimizer::fit(const char* options) } +bool RooMinimizer::fitFcn() const { + bool ret; + + switch (_fcnMode) { + case FcnMode::classic: { + ret = _theFitter->FitFCN(*static_cast(_fcn)); + break; + } + case FcnMode::gradient: { + auto thing = dynamic_cast(_fcn); + if (thing != nullptr) { + ret = _theFitter->FitFCN(*thing); + } else { + throw std::logic_error("In RooMinimizer::fitFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); + } + break; + } + } + + return ret; +} //////////////////////////////////////////////////////////////////////////////// @@ -317,7 +300,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 +329,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 +492,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 +521,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 +550,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,8 +799,61 @@ 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(static_cast(_fcn)); + } + case FcnMode::gradient: { + auto thing = dynamic_cast(_fcn); + if (thing != nullptr) { + return static_cast(thing); + } else { + throw std::logic_error("In RooMinimizer::fitterFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); + } + } + } + } +} +const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const +{ + if (getFitterMultiGenFcn()) { + switch (_fcnMode) { + case FcnMode::classic: { + return static_cast(static_cast(getFitterMultiGenFcn())); + } + case FcnMode::gradient: { + auto thing = dynamic_cast(getFitterMultiGenFcn()); + if (thing != nullptr) { + return static_cast(thing); + } else { + throw std::logic_error("In RooMinimizer::fitterFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); + } + } + } + } 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() ); +} + //////////////////////////////////////////////////////////////////////////////// /// Apply results of given external covariance matrix. i.e. propagate its errors diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 9ec4dd851a82a..220cb8e09ae40 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(testRooGradMinimizer RooGradMinimizer.cxx LIBRARIES RooFitCore) diff --git a/roofit/roofitcore/test/RooGradMinimizer.cxx b/roofit/roofitcore/test/RooGradMinimizer.cxx new file mode 100644 index 0000000000000..de3576a8074f2 --- /dev/null +++ b/roofit/roofitcore/test/RooGradMinimizer.cxx @@ -0,0 +1,847 @@ +/***************************************************************************** + * 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); + 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); + m1->setMinimizerType("Minuit2"); + + m1->setStrategy(0); + m1->setPrintLevel(100); + m1->setVerbose(kTRUE); + + m1->migrad(); +} + +/* +TEST(GradMinimizer, Gaussian2DVarToConst) { + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w = RooWorkspace(); + + w.factory("Gaussian::g1(x[-5,5],mu1[0,-3,3],sigma1[1])"); + RooRealVar *mu1 = w.var("mu1"); + + w.factory("Gaussian::g2(x[-5,5],mu2[4,-3,12],sigma2[2.5])"); + RooRealVar *mu2 = w.var("mu2"); + +// RooArgSet pdf_set = w.allPdfs(); + + // create event counts for all pdfs +// RooArgSet count_set; + + // ... for the gaussians + RooRealVar N1("Nsig1", "#signal events component 1", 5000, 0., 20000); + w.import(N1); + RooRealVar N2("Nsig2", "#signal events component 2", 5000, 0., 20000); + w.import(N2); +// count_set.add(*w.arg("Nsig1")); +// count_set.add(*w.arg("Nsig2")); + + RooAddPdf sum("sum", "2 gaussians", w.allPdfs(), RooArgSet(N1, N2)); + + auto x = w.var("x"); + RooDataSet *data = sum.generate(RooArgSet(*x), 10000); + mu1->setVal(-2.9); + mu2->setVal(1); + + auto nll = sum.createNLL(*data); + + // save initial values for the start of all minimizations + RooArgSet values = RooArgSet(RooArgSet(*mu1, *mu2, N1, N2, sum, *nll), w.allPdfs(), "all values"); + + 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(); + + values = *savedValues; + mu1->setConstant(kTRUE); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mu1_0 = mu1->getVal(); + double muerr1_0 = mu1->getError(); + double mu2_0 = mu2->getVal(); + double muerr2_0 = mu2->getError(); + + values = *savedValues; + mu1->setConstant(kFALSE); + + RooMinimizer m1(*nll); + m1.setMinimizerType("Minuit2"); + + m1.setStrategy(0); + m1.setPrintLevel(-1); + + m1.migrad(); + + values = *savedValues; + mu1->setConstant(kTRUE); + + m1.migrad(); + + RooFitResult *m1result = m1.lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mu1_1 = mu1->getVal(); + double muerr1_1 = mu1->getError(); + double mu2_1 = mu2->getVal(); + double muerr2_1 = mu2->getError(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(mu1_0, mu1_1); + EXPECT_EQ(muerr1_0, muerr1_1); + EXPECT_EQ(mu2_0, mu2_1); + EXPECT_EQ(muerr2_0, muerr2_1); + EXPECT_EQ(edm0, edm1); +} + +TEST(GradMinimizer, Gaussian2DConstToVar) { + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + // produce the same random stuff every time + RooRandom::randomGenerator()->SetSeed(1); + + RooWorkspace w = RooWorkspace(); + + w.factory("Gaussian::g1(x[-5,5],mu1[0,-3,3],sigma1[1])"); + RooRealVar *mu1 = w.var("mu1"); + + w.factory("Gaussian::g2(x[-5,5],mu2[4,-3,12],sigma2[2.5])"); + RooRealVar *mu2 = w.var("mu2"); + + RooArgSet pdf_set = w.allPdfs(); + + // create event counts for all pdfs + RooArgSet count_set; + + // ... for the gaussians + RooRealVar N1("Nsig1", "#signal events component 1", 5000, 0., 20000); + w.import(N1); + RooRealVar N2("Nsig2", "#signal events component 2", 5000, 0., 20000); + w.import(N2); + count_set.add(*w.arg("Nsig1")); + count_set.add(*w.arg("Nsig2")); + + RooAddPdf sum("sum", "2 gaussians", pdf_set, count_set); + + auto x = w.var("x"); + RooDataSet *data = sum.generate(RooArgSet(*x), 10000); + mu1->setVal(-2.9); + mu2->setVal(1); + + auto nll = sum.createNLL(*data); + + // save initial values for the start of all minimizations + RooArgSet values = RooArgSet(*mu1, *mu2, sum, *nll); + + 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); + + mu1->setConstant(kTRUE); + + m0.migrad(); + + values = *savedValues; + mu1->setConstant(kFALSE); + + m0.migrad(); + + RooFitResult *m0result = m0.lastMinuitFit(); + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mu1_0 = mu1->getVal(); + double muerr1_0 = mu1->getError(); + double mu2_0 = mu2->getVal(); + double muerr2_0 = mu2->getError(); + + values = *savedValues; + mu1->setConstant(kTRUE); + + RooMinimizer m1(*nll); + m1.setMinimizerType("Minuit2"); + + m1.setStrategy(0); + m1.setPrintLevel(-1); + + m1.migrad(); + + values = *savedValues; + mu1->setConstant(kFALSE); + + m1.migrad(); + + RooFitResult *m1result = m1.lastMinuitFit(); + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mu1_1 = mu1->getVal(); + double muerr1_1 = mu1->getError(); + double mu2_1 = mu2->getVal(); + double muerr2_1 = mu2->getError(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(mu1_0, mu1_1); + EXPECT_EQ(muerr1_0, muerr1_1); + EXPECT_EQ(mu2_0, mu2_1); + EXPECT_EQ(muerr2_0, muerr2_1); + EXPECT_EQ(edm0, edm1); +} +*/ + +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())); + + 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); + + 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); + + 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); + 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); + 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..3bafe17eff717 --- /dev/null +++ b/roofit/roofitcore/test/test_lib.h @@ -0,0 +1,185 @@ +/***************************************************************************** + * 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 + +#include // make_unique + +#include "RooWorkspace.h" +#include "RooRandom.h" +#include "RooAddPdf.h" +#include "RooDataSet.h" +#include "RooRealVar.h" // for the dynamic cast to have a complete type + + +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)); +} + + +class Hex { +public: + explicit Hex(double n) : number_(n) {} + operator double() const { return number_; } + bool operator==(const Hex& other) { + return double(*this) == double(other); + } + +private: + double number_; +}; + +::std::ostream& operator<<(::std::ostream& os, const Hex& hex) { + return os << std::hexfloat << double(hex) << std::defaultfloat; // whatever needed to print bar to os +} + + +#endif //ROOT_TEST_LIB_H From 885d6ce3c7a28f190f8aed23644bf4cfa11f327f Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Thu, 15 Jul 2021 17:02:43 +0200 Subject: [PATCH 2/7] apply review comments From comments by @hageboeck on PR #8596 (https://github.com/root-project/root/pull/8596): - Include reordering - Camel case member function names - Inline simple functions --- roofit/roofitcore/inc/RooGradMinimizerFcn.h | 62 +++++------ roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 102 ++++-------------- 2 files changed, 51 insertions(+), 113 deletions(-) diff --git a/roofit/roofitcore/inc/RooGradMinimizerFcn.h b/roofit/roofitcore/inc/RooGradMinimizerFcn.h index 87605092de2fd..e2619f0f943f8 100644 --- a/roofit/roofitcore/inc/RooGradMinimizerFcn.h +++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h @@ -16,58 +16,55 @@ #ifndef ROO_GRAD_MINIMIZER_FCN #define ROO_GRAD_MINIMIZER_FCN -#include +#include "RooAbsMinimizerFcn.h" #include "Minuit2/MnStrategy.h" -#include "Minuit2/MnMatrix.h" // MnAlgebraicVector -#include "Math/IFunction.h" // IMultiGradFunction +#include "Minuit2/MnMatrix.h" // MnAlgebraicVector +#include "Math/IFunction.h" // IMultiGradFunction #include "Fit/ParameterSettings.h" #include "Minuit2/NumericalDerivator.h" -#include "RooAbsMinimizerFcn.h" +#include class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbsMinimizerFcn { public: RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose = false); - RooGradMinimizerFcn(const RooGradMinimizerFcn& other); + RooGradMinimizerFcn(const RooGradMinimizerFcn &other); ROOT::Math::IMultiGradFunction *Clone() const override; - ROOT::Minuit2::MnStrategy get_strategy() const; - double get_error_def() const; - void set_strategy(int istrat); + 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 synchronize_gradient_parameter_settings(std::vector ¶meter_settings) const; - - bool returnsInMinuit2ParameterSpace() const override; - - unsigned int NDim() const override; + void synchronizeGradientParameterSettings(std::vector ¶meter_settings) const; - void set_step_tolerance(double step_tolerance) const; - void set_grad_tolerance(double grad_tolerance) const; - void set_ncycles(unsigned int ncycles) const; - void set_error_level(double error_level) const; - - std::string getFunctionName() const override; - std::string getFunctionTitle() const override; - - void setOffsetting(Bool_t flag) override; + 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.set_error_level(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 run_derivator(unsigned int i_component) const; + void runDerivator(unsigned int i_component) const; - void reset_has_been_calculated_flags() const; - bool sync_parameter(double x, std::size_t ix) const; - bool sync_parameters(const double *x) const; + void resetHasBeenCalculatedFlags() const; + bool syncParameter(double x, std::size_t ix) const; + bool syncParameters(const double *x) const; - void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override; + inline void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override + { + _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); + } public: - enum class GradientCalculatorMode { - ExactlyMinuit2, AlmostMinuit2 - }; + enum class GradientCalculatorMode { ExactlyMinuit2, AlmostMinuit2 }; private: // IMultiGradFunction overrides @@ -79,14 +76,11 @@ class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbs 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; - -public: - // for debugging, wraps ROOT::Math::IMultiGradFunction::Gradient, can be used for further inspection: - void Gradient(const double *x, double *grad) const override; }; #endif diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx index 08c259d186810..539c1d021dc85 100644 --- a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -21,25 +21,21 @@ #include "RooGradMinimizerFcn.h" -#include - -#include "Riostream.h" - -#include "TIterator.h" -#include "TClass.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(*funct->getParameters(RooArgSet())), context, verbose), @@ -49,9 +45,9 @@ RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *contex // TODO: added "parameters" after rewrite in april 2020, check if correct auto parameters = _context->fitter()->Config().ParamsSettings(); synchronizeParameterSettings(parameters, kTRUE, verbose); - synchronize_gradient_parameter_settings(parameters); - set_strategy(ROOT::Math::MinimizerOptions::DefaultStrategy()); - set_error_level(ROOT::Math::MinimizerOptions::DefaultErrorDef()); + synchronizeGradientParameterSettings(parameters); + setStrategy(ROOT::Math::MinimizerOptions::DefaultStrategy()); + setErrorLevel(ROOT::Math::MinimizerOptions::DefaultErrorDef()); } RooGradMinimizerFcn::RooGradMinimizerFcn(const RooGradMinimizerFcn &other) @@ -65,7 +61,7 @@ ROOT::Math::IMultiGradFunction *RooGradMinimizerFcn::Clone() const return new RooGradMinimizerFcn(*this); } -void RooGradMinimizerFcn::synchronize_gradient_parameter_settings( +void RooGradMinimizerFcn::synchronizeGradientParameterSettings( std::vector ¶meter_settings) const { _gradf.SetInitialGradient(_context->getMultiGenFcn(), parameter_settings, _grad); @@ -149,7 +145,7 @@ double RooGradMinimizerFcn::DoEval(const double *x) const return fvalue; } -void RooGradMinimizerFcn::reset_has_been_calculated_flags() const +void RooGradMinimizerFcn::resetHasBeenCalculatedFlags() const { for (auto it = has_been_calculated.begin(); it != has_been_calculated.end(); ++it) { *it = false; @@ -157,7 +153,7 @@ void RooGradMinimizerFcn::reset_has_been_calculated_flags() const none_have_been_calculated = true; } -bool RooGradMinimizerFcn::sync_parameter(double x, std::size_t ix) const +bool RooGradMinimizerFcn::syncParameter(double x, std::size_t ix) const { bool parameter_has_changed = (_grad_params[ix] != x); @@ -168,14 +164,14 @@ bool RooGradMinimizerFcn::sync_parameter(double x, std::size_t ix) const SetPdfParamVal(ix, x); if (!none_have_been_calculated) { - reset_has_been_calculated_flags(); + resetHasBeenCalculatedFlags(); } } return parameter_has_changed; } -bool RooGradMinimizerFcn::sync_parameters(const double *x) const +bool RooGradMinimizerFcn::syncParameters(const double *x) const { bool has_been_synced = false; @@ -193,13 +189,13 @@ bool RooGradMinimizerFcn::sync_parameters(const double *x) const } if (has_been_synced) { - reset_has_been_calculated_flags(); + resetHasBeenCalculatedFlags(); } return has_been_synced; } -void RooGradMinimizerFcn::run_derivator(unsigned int i_component) const +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]) { @@ -214,81 +210,29 @@ void RooGradMinimizerFcn::run_derivator(unsigned int i_component) const double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int i_component) const { - sync_parameters(x); - run_derivator(i_component); + syncParameters(x); + runDerivator(i_component); return _grad[i_component].derivative; } -bool RooGradMinimizerFcn::returnsInMinuit2ParameterSpace() const -{ - return true; -} - -unsigned int RooGradMinimizerFcn::NDim() const -{ - return getNDim(); -} - //////////////////////////////////////////////////////////////////////////////// -void RooGradMinimizerFcn::set_strategy(int istrat) +void RooGradMinimizerFcn::setStrategy(int istrat) { assert(istrat >= 0); ROOT::Minuit2::MnStrategy strategy(static_cast(istrat)); - set_step_tolerance(strategy.GradientStepTolerance()); - set_grad_tolerance(strategy.GradientTolerance()); - set_ncycles(strategy.GradientNCycles()); + 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); - synchronize_gradient_parameter_settings(parameters); - set_strategy(_context->fitter()->Config().MinimizerOptions().Strategy()); - set_error_level(_context->fitter()->Config().MinimizerOptions().ErrorDef()); + synchronizeGradientParameterSettings(parameters); + setStrategy(_context->fitter()->Config().MinimizerOptions().Strategy()); + setErrorLevel(_context->fitter()->Config().MinimizerOptions().ErrorDef()); return returnee; } - -void RooGradMinimizerFcn::setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) -{ - _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); -} - -void RooGradMinimizerFcn::set_step_tolerance(double step_tolerance) const -{ - _gradf.set_step_tolerance(step_tolerance); -} -void RooGradMinimizerFcn::set_grad_tolerance(double grad_tolerance) const -{ - _gradf.set_grad_tolerance(grad_tolerance); -} -void RooGradMinimizerFcn::set_ncycles(unsigned int ncycles) const -{ - _gradf.set_ncycles(ncycles); -} -void RooGradMinimizerFcn::set_error_level(double error_level) const -{ - _gradf.set_error_level(error_level); -} - -std::string RooGradMinimizerFcn::getFunctionName() const -{ - return _funct->GetName(); -} - -std::string RooGradMinimizerFcn::getFunctionTitle() const -{ - return _funct->GetTitle(); -} - -void RooGradMinimizerFcn::setOffsetting(Bool_t flag) -{ - _funct->enableOffsetting(flag); -} - -void RooGradMinimizerFcn::Gradient(const double *x, double *grad) const -{ - ROOT::Math::IMultiGradFunction::Gradient(x, grad); -} From cab5a8fa33e1207a9ce1e094500f5cfa61ddeb73 Mon Sep 17 00:00:00 2001 From: Patrick Bos Date: Thu, 15 Jul 2021 17:51:20 +0200 Subject: [PATCH 3/7] Apply suggestions from code review Co-authored-by: Jonas Rembser Co-authored-by: Stephan Hageboeck --- roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx index 539c1d021dc85..a5516620ed0f3 100644 --- a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -38,7 +38,7 @@ #include RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose) - : RooAbsMinimizerFcn(RooArgList(*funct->getParameters(RooArgSet())), context, verbose), + : RooAbsMinimizerFcn(RooArgList( * std::unique_ptr(funct->getParameters({})) ), context, verbose), _grad(getNDim()), _grad_params(getNDim()), _funct(funct), has_been_calculated(getNDim()) { From 8002b33a0215d96a2f11a090638baf05f324e5a4 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Thu, 15 Jul 2021 20:37:54 +0200 Subject: [PATCH 4/7] apply review comments From comments by @guitargeek and @hageboeck on PR #8596 (https://github.com/root-project/root/pull/8596): In RooMinimizer: - Removed template parameter from RooMinimizer ctor. - Remove dynamic_cast nullptr checks when FcnMode is already known. - Also replace static with dynamic_casts when downcasting; fixes clang tidy warning. - Fix compile warnings due to not returning anything at end of non-void function with switch statements. Other places: - In test_lib: remove Hex number printing class and reorder includes. - Renamed test to testRooGradMinimizerFcn(.cxx) and cleaned up old commented out code. --- roofit/roofitcore/inc/RooMinimizer.h | 11 +- roofit/roofitcore/src/RooMinimizer.cxx | 99 ++++++-- roofit/roofitcore/test/CMakeLists.txt | 2 +- ...imizer.cxx => testRooGradMinimizerFcn.cxx} | 216 +----------------- roofit/roofitcore/test/test_lib.h | 23 +- 5 files changed, 87 insertions(+), 264 deletions(-) rename roofit/roofitcore/test/{RooGradMinimizer.cxx => testRooGradMinimizerFcn.cxx} (77%) diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 6eb930aa733a7..a281395e74a7f 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -50,10 +50,9 @@ class RooMinimizer : public TObject { public: enum class FcnMode { classic, gradient }; - RooMinimizer(RooAbsReal& function) ; - template - static std::unique_ptr create(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 } ; @@ -127,9 +126,6 @@ class RooMinimizer : public TObject { bool fitFcn() const; private: - template - RooMinimizer(RooAbsReal &function, MinimizerFcn* /* used only for template deduction */); - Int_t _printLevel = 1; Int_t _status = -99; Bool_t _profile = kFALSE; @@ -154,5 +150,4 @@ class RooMinimizer : public TObject { ClassDef(RooMinimizer,0) // RooFit interface to ROOT::Fit::Fitter } ; - #endif diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index df1bd6d8f1316..4b7c9c285e0a9 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -95,11 +95,66 @@ void RooMinimizer::cleanup() //////////////////////////////////////////////////////////////////////////////// -/// Constructor - -RooMinimizer::RooMinimizer(RooAbsReal &function) : - RooMinimizer::RooMinimizer(function, static_cast(nullptr)) +/// 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. + +RooMinimizer::RooMinimizer(RooAbsReal &function, FcnMode fcnMode) : _fcnMode(fcnMode) { + 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); + 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); } @@ -265,18 +320,16 @@ bool RooMinimizer::fitFcn() const { switch (_fcnMode) { case FcnMode::classic: { - ret = _theFitter->FitFCN(*static_cast(_fcn)); + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); break; } case FcnMode::gradient: { - auto thing = dynamic_cast(_fcn); - if (thing != nullptr) { - ret = _theFitter->FitFCN(*thing); - } else { - throw std::logic_error("In RooMinimizer::fitFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); - } + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); break; } + default: { + throw std::logic_error("In RooMinimizer::fitFcn: _fcnMode has an unsupported value!"); + } } return ret; @@ -812,15 +865,13 @@ ROOT::Math::IMultiGenFunction* RooMinimizer::getMultiGenFcn() const } else { switch (_fcnMode) { case FcnMode::classic: { - return static_cast(static_cast(_fcn)); + return static_cast(dynamic_cast(_fcn)); } case FcnMode::gradient: { - auto thing = dynamic_cast(_fcn); - if (thing != nullptr) { - return static_cast(thing); - } else { - throw std::logic_error("In RooMinimizer::fitterFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); - } + return static_cast(dynamic_cast(_fcn)); + } + default: { + throw std::logic_error("In RooMinimizer::getMultiGenFcn: _fcnMode has an unsupported value!"); } } } @@ -832,15 +883,13 @@ const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const if (getFitterMultiGenFcn()) { switch (_fcnMode) { case FcnMode::classic: { - return static_cast(static_cast(getFitterMultiGenFcn())); + return static_cast(dynamic_cast(getFitterMultiGenFcn())); } case FcnMode::gradient: { - auto thing = dynamic_cast(getFitterMultiGenFcn()); - if (thing != nullptr) { - return static_cast(thing); - } else { - throw std::logic_error("In RooMinimizer::fitterFcn: Minimizer fcnMode was set to gradient, but fit function type does not match RooGradMinimizerFcn!"); - } + return static_cast(dynamic_cast(getFitterMultiGenFcn())); + } + default: { + throw std::logic_error("In RooMinimizer::fitterFcn: _fcnMode has an unsupported value!"); } } } else { diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 220cb8e09ae40..07ae8638aba90 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -40,5 +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(testRooGradMinimizer RooGradMinimizer.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore) diff --git a/roofit/roofitcore/test/RooGradMinimizer.cxx b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx similarity index 77% rename from roofit/roofitcore/test/RooGradMinimizer.cxx rename to roofit/roofitcore/test/testRooGradMinimizerFcn.cxx index de3576a8074f2..2ce38eb71a800 100644 --- a/roofit/roofitcore/test/RooGradMinimizer.cxx +++ b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx @@ -71,7 +71,7 @@ TEST(GradMinimizer, Gaussian1D) *values = *savedValues; - std::unique_ptr m1 = RooMinimizer::create(*nll); + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); m1->setMinimizerType("Minuit2"); m1->setStrategy(0); @@ -132,7 +132,7 @@ TEST(GradMinimizerDebugging, DISABLED_Gaussian1DGradMinimizer) // when c++17 support arrives, change to this: // auto [nll, _] = generate_1D_gaussian_pdf_nll(w, 10000); - std::unique_ptr m1 = RooMinimizer::create(*nll); + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); m1->setMinimizerType("Minuit2"); m1->setStrategy(0); @@ -142,208 +142,6 @@ TEST(GradMinimizerDebugging, DISABLED_Gaussian1DGradMinimizer) m1->migrad(); } -/* -TEST(GradMinimizer, Gaussian2DVarToConst) { - RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); - - // produce the same random stuff every time - RooRandom::randomGenerator()->SetSeed(1); - - RooWorkspace w = RooWorkspace(); - - w.factory("Gaussian::g1(x[-5,5],mu1[0,-3,3],sigma1[1])"); - RooRealVar *mu1 = w.var("mu1"); - - w.factory("Gaussian::g2(x[-5,5],mu2[4,-3,12],sigma2[2.5])"); - RooRealVar *mu2 = w.var("mu2"); - -// RooArgSet pdf_set = w.allPdfs(); - - // create event counts for all pdfs -// RooArgSet count_set; - - // ... for the gaussians - RooRealVar N1("Nsig1", "#signal events component 1", 5000, 0., 20000); - w.import(N1); - RooRealVar N2("Nsig2", "#signal events component 2", 5000, 0., 20000); - w.import(N2); -// count_set.add(*w.arg("Nsig1")); -// count_set.add(*w.arg("Nsig2")); - - RooAddPdf sum("sum", "2 gaussians", w.allPdfs(), RooArgSet(N1, N2)); - - auto x = w.var("x"); - RooDataSet *data = sum.generate(RooArgSet(*x), 10000); - mu1->setVal(-2.9); - mu2->setVal(1); - - auto nll = sum.createNLL(*data); - - // save initial values for the start of all minimizations - RooArgSet values = RooArgSet(RooArgSet(*mu1, *mu2, N1, N2, sum, *nll), w.allPdfs(), "all values"); - - 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(); - - values = *savedValues; - mu1->setConstant(kTRUE); - - m0.migrad(); - - RooFitResult *m0result = m0.lastMinuitFit(); - double minNll0 = m0result->minNll(); - double edm0 = m0result->edm(); - double mu1_0 = mu1->getVal(); - double muerr1_0 = mu1->getError(); - double mu2_0 = mu2->getVal(); - double muerr2_0 = mu2->getError(); - - values = *savedValues; - mu1->setConstant(kFALSE); - - RooMinimizer m1(*nll); - m1.setMinimizerType("Minuit2"); - - m1.setStrategy(0); - m1.setPrintLevel(-1); - - m1.migrad(); - - values = *savedValues; - mu1->setConstant(kTRUE); - - m1.migrad(); - - RooFitResult *m1result = m1.lastMinuitFit(); - double minNll1 = m1result->minNll(); - double edm1 = m1result->edm(); - double mu1_1 = mu1->getVal(); - double muerr1_1 = mu1->getError(); - double mu2_1 = mu2->getVal(); - double muerr2_1 = mu2->getError(); - - EXPECT_EQ(minNll0, minNll1); - EXPECT_EQ(mu1_0, mu1_1); - EXPECT_EQ(muerr1_0, muerr1_1); - EXPECT_EQ(mu2_0, mu2_1); - EXPECT_EQ(muerr2_0, muerr2_1); - EXPECT_EQ(edm0, edm1); -} - -TEST(GradMinimizer, Gaussian2DConstToVar) { - RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); - - // produce the same random stuff every time - RooRandom::randomGenerator()->SetSeed(1); - - RooWorkspace w = RooWorkspace(); - - w.factory("Gaussian::g1(x[-5,5],mu1[0,-3,3],sigma1[1])"); - RooRealVar *mu1 = w.var("mu1"); - - w.factory("Gaussian::g2(x[-5,5],mu2[4,-3,12],sigma2[2.5])"); - RooRealVar *mu2 = w.var("mu2"); - - RooArgSet pdf_set = w.allPdfs(); - - // create event counts for all pdfs - RooArgSet count_set; - - // ... for the gaussians - RooRealVar N1("Nsig1", "#signal events component 1", 5000, 0., 20000); - w.import(N1); - RooRealVar N2("Nsig2", "#signal events component 2", 5000, 0., 20000); - w.import(N2); - count_set.add(*w.arg("Nsig1")); - count_set.add(*w.arg("Nsig2")); - - RooAddPdf sum("sum", "2 gaussians", pdf_set, count_set); - - auto x = w.var("x"); - RooDataSet *data = sum.generate(RooArgSet(*x), 10000); - mu1->setVal(-2.9); - mu2->setVal(1); - - auto nll = sum.createNLL(*data); - - // save initial values for the start of all minimizations - RooArgSet values = RooArgSet(*mu1, *mu2, sum, *nll); - - 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); - - mu1->setConstant(kTRUE); - - m0.migrad(); - - values = *savedValues; - mu1->setConstant(kFALSE); - - m0.migrad(); - - RooFitResult *m0result = m0.lastMinuitFit(); - double minNll0 = m0result->minNll(); - double edm0 = m0result->edm(); - double mu1_0 = mu1->getVal(); - double muerr1_0 = mu1->getError(); - double mu2_0 = mu2->getVal(); - double muerr2_0 = mu2->getError(); - - values = *savedValues; - mu1->setConstant(kTRUE); - - RooMinimizer m1(*nll); - m1.setMinimizerType("Minuit2"); - - m1.setStrategy(0); - m1.setPrintLevel(-1); - - m1.migrad(); - - values = *savedValues; - mu1->setConstant(kFALSE); - - m1.migrad(); - - RooFitResult *m1result = m1.lastMinuitFit(); - double minNll1 = m1result->minNll(); - double edm1 = m1result->edm(); - double mu1_1 = mu1->getVal(); - double muerr1_1 = mu1->getError(); - double mu2_1 = mu2->getVal(); - double muerr2_1 = mu2->getError(); - - EXPECT_EQ(minNll0, minNll1); - EXPECT_EQ(mu1_0, mu1_1); - EXPECT_EQ(muerr1_0, muerr1_1); - EXPECT_EQ(mu2_0, mu2_1); - EXPECT_EQ(muerr2_0, muerr2_1); - EXPECT_EQ(edm0, edm1); -} -*/ - TEST(GradMinimizer, GaussianND) { // test RooMinimizer class with simple N-dimensional pdf @@ -405,7 +203,7 @@ TEST(GradMinimizer, GaussianND) // -------- - std::unique_ptr m1 = RooMinimizer::create(*(nll.get())); + std::unique_ptr m1 = RooMinimizer::create(*(nll.get()), RooMinimizer::FcnMode::gradient); m1->setStrategy(0); m1->setPrintLevel(-1); @@ -468,7 +266,7 @@ TEST(GradMinimizerReverse, GaussianND) // -------- - std::unique_ptr m0 = RooMinimizer::create(*nll); + std::unique_ptr m0 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); m0->setMinimizerType("Minuit2"); @@ -652,7 +450,7 @@ TEST(GradMinimizer, BranchingPDF) // -------- - std::unique_ptr m1 = RooMinimizer::create(*nll); + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); m1->setStrategy(0); m1->setPrintLevel(-1); @@ -769,7 +567,7 @@ TEST(GradMinimizerDebugging, DISABLED_BranchingPDFLoadFromWorkspace) all_values.Print("v"); - std::unique_ptr m1 = RooMinimizer::create(*nll); + std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); m1->setMinimizerType("Minuit2"); m1->setStrategy(0); @@ -840,7 +638,7 @@ TEST(GradMinimizerDebugging, DISABLED_BranchingPDFLoadFromWorkspaceGradMinimizer RooDataSet *data = static_cast(w.data("")); auto nll = sum.createNLL(*data); - std::unique_ptr m0 = RooMinimizer::create(*nll); + 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 index 3bafe17eff717..dae7ef3b50dd7 100644 --- a/roofit/roofitcore/test/test_lib.h +++ b/roofit/roofitcore/test/test_lib.h @@ -13,16 +13,14 @@ #ifndef ROOT_TEST_LIB_H #define ROOT_TEST_LIB_H -#include - -#include // make_unique - #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) { @@ -165,21 +163,4 @@ generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_ev } -class Hex { -public: - explicit Hex(double n) : number_(n) {} - operator double() const { return number_; } - bool operator==(const Hex& other) { - return double(*this) == double(other); - } - -private: - double number_; -}; - -::std::ostream& operator<<(::std::ostream& os, const Hex& hex) { - return os << std::hexfloat << double(hex) << std::defaultfloat; // whatever needed to print bar to os -} - - #endif //ROOT_TEST_LIB_H From 329903b71c335ab677337adcf7d4784983104330 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 16 Jul 2021 11:44:15 +0200 Subject: [PATCH 5/7] temporarily disable testRooGradMinimizerFcn In an upcoming PR, we will introduce changes to Minuit2 and mathcore, necessary for fixing these tests. --- roofit/roofitcore/test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 07ae8638aba90..9a1a54c81d4d1 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -40,5 +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) +#ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore) From bebd62c514f881ee7a66d70f1643be1cc0694dbc Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 16 Jul 2021 13:48:26 +0200 Subject: [PATCH 6/7] fix compilation errors after rebase --- roofit/roofitcore/inc/RooGradMinimizerFcn.h | 2 +- roofit/roofitcore/inc/RooMinimizer.h | 2 +- roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/roofit/roofitcore/inc/RooGradMinimizerFcn.h b/roofit/roofitcore/inc/RooGradMinimizerFcn.h index e2619f0f943f8..829fc228d45d5 100644 --- a/roofit/roofitcore/inc/RooGradMinimizerFcn.h +++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h @@ -46,7 +46,7 @@ class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbs 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.set_error_level(error_level); } + 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); } diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index a281395e74a7f..c42975663c8ed 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -147,7 +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 index a5516620ed0f3..55e5d66ab5010 100644 --- a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx @@ -38,7 +38,7 @@ #include RooGradMinimizerFcn::RooGradMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose) - : RooAbsMinimizerFcn(RooArgList( * std::unique_ptr(funct->getParameters({})) ), context, verbose), + : RooAbsMinimizerFcn(RooArgList( * std::unique_ptr(funct->getParameters(RooArgSet{})) ), context, verbose), _grad(getNDim()), _grad_params(getNDim()), _funct(funct), has_been_calculated(getNDim()) { @@ -201,8 +201,8 @@ void RooGradMinimizerFcn::runDerivator(unsigned int i_component) const if (!has_been_calculated[i_component]) { // Calculate the derivative etc for these parameters _grad[i_component] = - _gradf.partial_derivative(_context->getMultiGenFcn(), _grad_params.data(), - _context->fitter()->Config().ParamsSettings(), i_component, _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; } From c7fafc0682d1e25bd9af3563dd0f534661befa09 Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Fri, 16 Jul 2021 15:28:57 +0200 Subject: [PATCH 7/7] revert default change Minuit -> Minuit2 --- roofit/roofitcore/inc/RooMinimizer.h | 2 +- roofit/roofitcore/src/RooMinimizer.cxx | 8 +++++--- roofit/roofitcore/test/testRooGradMinimizerFcn.cxx | 2 ++ 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index c42975663c8ed..835e71b0aa791 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -138,7 +138,7 @@ class RooMinimizer : public TObject { TMatrixDSym *_extV = 0; RooAbsMinimizerFcn *_fcn; - std::string _minimizerType = "Minuit2"; + std::string _minimizerType = "Minuit"; FcnMode _fcnMode; static ROOT::Fit::Fitter *_theFitter ; diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 4b7c9c285e0a9..7bcdab4c708c7 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -21,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 MINUIT2. +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 @@ -124,6 +125,7 @@ RooMinimizer::RooMinimizer(RooAbsReal &function, FcnMode fcnMode) : _fcnMode(fcn } case FcnMode::gradient: { _fcn = new RooGradMinimizerFcn(&function, this, _verbose); + setMinimizerType("Minuit2"); break; } default: { @@ -254,8 +256,8 @@ class RooGradMinimizerFcn; void RooMinimizer::setMinimizerType(const char* type) { - if (dynamic_cast(_fcn) && strcmp(type, "Minuit2") != 0) { - throw std::invalid_argument("In RooMinimizer::setMinimizerType: only Minuit2 is supported when using RooGradMinimizerFcn!"); + 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; } diff --git a/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx index 2ce38eb71a800..3d1a8181612c8 100644 --- a/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx +++ b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx @@ -204,6 +204,7 @@ TEST(GradMinimizer, GaussianND) // -------- std::unique_ptr m1 = RooMinimizer::create(*(nll.get()), RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); m1->setStrategy(0); m1->setPrintLevel(-1); @@ -451,6 +452,7 @@ TEST(GradMinimizer, BranchingPDF) // -------- std::unique_ptr m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient); + m1->setMinimizerType("Minuit2"); m1->setStrategy(0); m1->setPrintLevel(-1);