diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 3ebb292c39441..bcba964cd9570 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -139,6 +139,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore RooMappedCategory.h RooMCIntegrator.h RooMCStudy.h + RooCladMinimizerFcn.h RooMinimizerFcn.h RooMinimizer.h RooMoment.h @@ -368,6 +369,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/RooMCIntegrator.cxx src/RooMCStudy.cxx src/RooMinimizer.cxx + src/RooCladMinimizerFcn.cxx src/RooMinimizerFcn.cxx src/RooMoment.cxx src/RooMPSentinel.cxx @@ -460,6 +462,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/TestStatistics/LikelihoodGradientWrapper.cxx src/TestStatistics/LikelihoodWrapper.cxx src/TestStatistics/LikelihoodSerial.cxx + src/TestStatistics/LikelihoodGradientClad.cxx src/TestStatistics/MinuitFcnGrad.cxx src/TestStatistics/RooAbsL.cxx src/TestStatistics/RooBinnedL.cxx diff --git a/roofit/roofitcore/inc/RooAbsReal.h b/roofit/roofitcore/inc/RooAbsReal.h index cbf53a3421370..11108e27b2e45 100644 --- a/roofit/roofitcore/inc/RooAbsReal.h +++ b/roofit/roofitcore/inc/RooAbsReal.h @@ -374,6 +374,8 @@ class RooAbsReal : public RooAbsArg { static void setHideOffset(Bool_t flag); static Bool_t hideOffset() ; + virtual void evaluateGradient(double*) const {} + protected: // Hook for objects with normalization-dependent parameters interpretation virtual void selectNormalization(const RooArgSet* depSet=0, Bool_t force=kFALSE) ; diff --git a/roofit/roofitcore/inc/RooCladMinimizerFcn.h b/roofit/roofitcore/inc/RooCladMinimizerFcn.h new file mode 100644 index 0000000000000..c1cd14b878903 --- /dev/null +++ b/roofit/roofitcore/inc/RooCladMinimizerFcn.h @@ -0,0 +1,71 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, Netherlands 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 roofit_roofitcore_RooCladMinimizerFcn_h +#define roofit_roofitcore_RooCladMinimizerFcn_h + +#include "Math/IFunction.h" +#include "Fit/ParameterSettings.h" +#include "Fit/FitResult.h" + +#include "RooAbsReal.h" +#include "RooArgList.h" + +#include +#include + +#include + +template class TMatrixTSym; +using TMatrixDSym = TMatrixTSym; + +// forward declaration +class RooMinimizer; + +// class RooCladMinimizerFcn : public RooAbsMinimizerFcn, public ROOT::Math::IBaseFunctionMultiDim { +class RooCladMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbsMinimizerFcn { + +public: + RooCladMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose = false); + RooCladMinimizerFcn(const RooCladMinimizerFcn &other); + ~RooCladMinimizerFcn() override; + + ROOT::Math::IBaseFunctionMultiDim *Clone() const override; + + /// IMultiGradFunction overrides necessary for Minuit + void Gradient(const double *x, double *grad) const override; + void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, + double *previous_gstep) const override; + + unsigned int NDim() const override { return getNDim(); } + + std::string getFunctionName() const override; + std::string getFunctionTitle() const override; + + void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override; + + void setOffsetting(Bool_t flag) override; + +private: + /// IMultiGradFunction overrides necessary for Minuit + double DoDerivative(const double *x, unsigned int icoord) const override; + double DoDerivativeWithPrevResult(const double *x, unsigned int i_component, double *previous_grad, + double *previous_g2, double *previous_gstep) const override; + + double DoEval(const double *x) const override; + + RooAbsReal *_funct; +}; + +#endif diff --git a/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h b/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h index ebc6e80b2f275..9117a469f1340 100644 --- a/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h +++ b/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h @@ -30,7 +30,7 @@ namespace TestStatistics { class RooAbsL; struct WrapperCalculationCleanFlags; -enum class LikelihoodGradientMode { multiprocess }; +enum class LikelihoodGradientMode { multiprocess, clad }; class LikelihoodGradientWrapper { public: diff --git a/roofit/roofitcore/inc/RooFit/TestStatistics/RooAbsL.h b/roofit/roofitcore/inc/RooFit/TestStatistics/RooAbsL.h index f521cd36a9135..81e15bced13c7 100644 --- a/roofit/roofitcore/inc/RooFit/TestStatistics/RooAbsL.h +++ b/roofit/roofitcore/inc/RooFit/TestStatistics/RooAbsL.h @@ -95,6 +95,8 @@ class RooAbsL { virtual ROOT::Math::KahanSum evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) = 0; + virtual void evaluateGradient(double*) const {} + // necessary from MinuitFcnGrad to reach likelihood properties: virtual RooArgSet *getParameters(); diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 2b8e44d9c5fca..d57e63d606ea8 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -48,7 +48,7 @@ class RooPlot ; class RooMinimizer : public TObject { public: - enum class FcnMode { classic, gradient, generic_wrapper }; + enum class FcnMode { classic, gradient, generic_wrapper, clad }; explicit RooMinimizer(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic); explicit RooMinimizer(std::shared_ptr likelihood, diff --git a/roofit/roofitcore/src/RooCladMinimizerFcn.cxx b/roofit/roofitcore/src/RooCladMinimizerFcn.cxx new file mode 100644 index 0000000000000..a1c0d8dd20ac3 --- /dev/null +++ b/roofit/roofitcore/src/RooCladMinimizerFcn.cxx @@ -0,0 +1,165 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * @(#)root/roofitcore:$Id$ + * Authors: * + * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it * + * PB, Patrick Bos, Netherlands 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) * + *****************************************************************************/ + +////////////////////////////////////////////////////////////////////////////// +/// \class RooCladMinimizerFcn +/// RooCladMinimizerFcn is an interface to the ROOT::Math::IBaseFunctionMultiDim, +/// a function that ROOT's minimisers use to carry out minimisations. +/// + +#include "RooCladMinimizerFcn.h" + +#include "RooAbsArg.h" +#include "RooAbsPdf.h" +#include "RooArgSet.h" +#include "RooRealVar.h" +#include "RooAbsRealLValue.h" +#include "RooMsgService.h" +#include "RooMinimizer.h" +#include "RooNaNPacker.h" + +#include "TClass.h" +#include "TMatrixDSym.h" + +#include +#include + +using namespace std; + + +namespace { + +// Helper function that wraps RooAbsArg::getParameters and directly returns the +// output RooArgSet. To be used in the initializer list of the RooCladMinimizerFcn +// constructor. +RooArgSet getParameters(RooAbsReal const& funct) { + RooArgSet out; + funct.getParameters(nullptr, out); + return out; +} + +} // namespace + + +RooCladMinimizerFcn::RooCladMinimizerFcn(RooAbsReal *funct, RooMinimizer* context, + bool verbose) : + RooAbsMinimizerFcn(getParameters(*funct), context, verbose), _funct(funct) +{} + + + +RooCladMinimizerFcn::RooCladMinimizerFcn(const RooCladMinimizerFcn& other) : ROOT::Math::IMultiGradFunction(other), RooAbsMinimizerFcn(other), + _funct(other._funct) +{} + + +RooCladMinimizerFcn::~RooCladMinimizerFcn() +{} + + +ROOT::Math::IBaseFunctionMultiDim* RooCladMinimizerFcn::Clone() const +{ + return new RooCladMinimizerFcn(*this) ; +} + +void RooCladMinimizerFcn::setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) +{ + _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); +} + +/// Evaluate function given the parameters in `x`. +double RooCladMinimizerFcn::DoEval(const double *x) const { + + // Set the parameter values for this iteration + for (unsigned index = 0; index < _nDim; index++) { + if (_logfile) (*_logfile) << x[index] << " " ; + SetPdfParamVal(index,x[index]); + } + + // Calculate the function for these parameters + RooAbsReal::setHideOffset(kFALSE) ; + double fvalue = _funct->getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + + if (!std::isfinite(fvalue) || RooAbsReal::numEvalErrors() > 0 || fvalue > 1e30) { + printEvalErrors(); + RooAbsReal::clearEvalErrorLog() ; + _numBadNLL++ ; + + if (_doEvalErrorWall) { + const double badness = RooNaNPacker::unpackNaN(fvalue); + fvalue = (std::isfinite(_maxFCN) ? _maxFCN : 0.) + _recoverFromNaNStrength * badness; + } + } else { + if (_evalCounter > 0 && _evalCounter == _numBadNLL) { + // This is the first time we get a valid function value; while before, the + // function was always invalid. For invalid cases, we returned values > 0. + // Now, we offset valid values such that they are < 0. + _funcOffset = -fvalue; + } + fvalue += _funcOffset; + _maxFCN = std::max(fvalue, _maxFCN); + } + + // Optional logging + if (_logfile) + (*_logfile) << setprecision(15) << fvalue << setprecision(4) << endl; + if (_verbose) { + cout << "\nprevFCN" << (_funct->isOffsetting()?"-offset":"") << " = " << setprecision(10) + << fvalue << setprecision(4) << " " ; + cout.flush() ; + } + + _evalCounter++ ; + + return fvalue; +} + +std::string RooCladMinimizerFcn::getFunctionName() const +{ + return _funct->GetName(); +} + +std::string RooCladMinimizerFcn::getFunctionTitle() const +{ + return _funct->GetTitle(); +} + +void RooCladMinimizerFcn::setOffsetting(Bool_t flag) +{ + _funct->enableOffsetting(flag); +} + +void RooCladMinimizerFcn::Gradient(const double *x, double *grad) const { + DoEval(x); + _funct->evaluateGradient(grad); +} + +void RooCladMinimizerFcn::GradientWithPrevResult(const double *x, double *grad, double * /*previous_grad*/, double * /*previous_g2*/, + double * /*previous_gstep*/) const { + Gradient(x, grad); +} + +double RooCladMinimizerFcn::DoDerivative(const double *x, unsigned int icoord) const { + double grad[_nDim]; + Gradient(x, grad); + return grad[icoord]; +} + +double RooCladMinimizerFcn::DoDerivativeWithPrevResult(const double *x, unsigned int i_component, double *previous_grad, + double *previous_g2, double *previous_gstep) const { + double grad[_nDim]; + Gradient(x, grad); + return grad[i_component]; +} diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index c81702588c174..0c31928bfb394 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -51,6 +51,7 @@ automatic PDF optimization. #include "RooMsgService.h" #include "RooPlot.h" #include "RooMinimizerFcn.h" +#include "RooCladMinimizerFcn.h" #include "RooGradMinimizerFcn.h" #include "RooFitResult.h" #include "TestStatistics/MinuitFcnGrad.h" @@ -108,6 +109,10 @@ RooMinimizer::RooMinimizer(RooAbsReal &function, FcnMode fcnMode) : _fcnMode(fcn _fcn = new RooMinimizerFcn(&function, this, _verbose); break; } + case FcnMode::clad: { + _fcn = new RooCladMinimizerFcn(&function, this, _verbose); + break; + } case FcnMode::gradient: { _fcn = new RooGradMinimizerFcn(&function, this, _verbose); setMinimizerType("Minuit2"); @@ -297,6 +302,10 @@ bool RooMinimizer::fitFcn() const { ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); break; } + case FcnMode::clad: { + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); + break; + } case FcnMode::gradient: { ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); break; @@ -835,6 +844,9 @@ ROOT::Math::IMultiGenFunction* RooMinimizer::getMultiGenFcn() const case FcnMode::classic: { return static_cast(dynamic_cast(_fcn)); } + case FcnMode::clad: { + return static_cast(dynamic_cast(_fcn)); + } case FcnMode::gradient: { return static_cast(dynamic_cast(_fcn)); } @@ -859,6 +871,9 @@ const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const case FcnMode::gradient: { return static_cast(dynamic_cast(getFitterMultiGenFcn())); } + case FcnMode::clad: { + return static_cast(dynamic_cast(getFitterMultiGenFcn())); + } case FcnMode::generic_wrapper: { return static_cast(dynamic_cast(getFitterMultiGenFcn())); } diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.cxx new file mode 100644 index 0000000000000..ad13ccc6540fc --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.cxx @@ -0,0 +1,87 @@ +/* + * Project: RooFit + * Authors: + * PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl + * + * Copyright (c) 2021, CERN + * + * 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 "LikelihoodGradientClad.h" + +#include "RooFit/MultiProcess/JobManager.h" +#include "RooFit/MultiProcess/Messenger.h" +#include "RooFit/MultiProcess/Queue.h" +#include "RooMsgService.h" +#include "RooMinimizer.h" + +#include "Minuit2/MnStrategy.h" + +namespace RooFit { +namespace TestStatistics { + +LikelihoodGradientClad::LikelihoodGradientClad(std::shared_ptr likelihood, + std::shared_ptr calculation_is_clean, + std::size_t N_dim, RooMinimizer *minimizer) + : LikelihoodGradientWrapper(std::move(likelihood), std::move(calculation_is_clean), N_dim, minimizer) +{ + // Note to future maintainers: take care when storing the minimizer_fcn pointer. The + // RooAbsMinimizerFcn subclasses may get cloned inside MINUIT, which means the pointer + // should also somehow be updated in this class. + minuit_internal_x_.reserve(N_dim); + + _grad.resize(N_dim); +} + +LikelihoodGradientClad::LikelihoodGradientClad(const LikelihoodGradientClad &other) + : LikelihoodGradientWrapper(other), + minuit_internal_x_(other.minuit_internal_x_) +{ + _grad.resize(minimizer_->getNPar()); +} + +LikelihoodGradientClad *LikelihoodGradientClad::clone() const +{ + return new LikelihoodGradientClad(*this); +} + +void LikelihoodGradientClad::synchronizeParameterSettings( + const std::vector ¶meter_settings) +{ + LikelihoodGradientWrapper::synchronizeParameterSettings(parameter_settings); +} + +void LikelihoodGradientClad::synchronizeParameterSettings( + ROOT::Math::IMultiGenFunction *function, const std::vector ¶meter_settings) +{ +} + +void LikelihoodGradientClad::fillGradient(double *grad) +{ + if (!calculation_is_clean_->gradient) { + likelihood_->evaluateGradient(_grad.data()); + } + + for(std::size_t i = 0; i < minimizer_->getNPar(); ++i) grad[i] = _grad[i]; +} + +void LikelihoodGradientClad::fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, + double *previous_gstep) +{ + if (!calculation_is_clean_->gradient) { + likelihood_->evaluateGradient(_grad.data()); + } + + for(std::size_t i = 0; i < minimizer_->getNPar(); ++i) grad[i] = _grad[i]; +} + +void LikelihoodGradientClad::updateMinuitInternalParameterValues(const std::vector &minuit_internal_x) +{ + minuit_internal_x_ = minuit_internal_x; +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.h b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.h new file mode 100644 index 0000000000000..5310049e20a18 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientClad.h @@ -0,0 +1,63 @@ +/* + * Project: RooFit + * Authors: + * PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl + * + * Copyright (c) 2021, CERN + * + * 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_ROOFIT_TESTSTATISTICS_LikelihoodGradientClad +#define ROOT_ROOFIT_TESTSTATISTICS_LikelihoodGradientClad + +#include "RooFit/MultiProcess/Job.h" +#include "RooFit/TestStatistics/LikelihoodGradientWrapper.h" + +#include "Math/MinimizerOptions.h" +#include "Minuit2/NumericalDerivator.h" +#include "Minuit2/MnMatrix.h" + +#include + +namespace RooFit { +namespace TestStatistics { + +class LikelihoodGradientClad : public LikelihoodGradientWrapper { +public: + LikelihoodGradientClad(std::shared_ptr likelihood, + std::shared_ptr calculation_is_clean, std::size_t N_dim, + RooMinimizer *minimizer); + LikelihoodGradientClad *clone() const override; + LikelihoodGradientClad(const LikelihoodGradientClad &other); + + void fillGradient(double *grad) override; + void fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, + double *previous_gstep) override; + + enum class GradientCalculatorMode { ExactlyMinuit2, AlmostMinuit2 }; + +private: + void synchronizeParameterSettings(ROOT::Math::IMultiGenFunction *function, + const std::vector ¶meter_settings) override; + // this overload must also be overridden here so that the one above doesn't trigger a overloaded-virtual warning: + void synchronizeParameterSettings(const std::vector ¶meter_settings) override; + + void synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions &options) override {} + + void updateMinuitInternalParameterValues(const std::vector &minuit_internal_x) override; + + bool usesMinuitInternalValues() override { return false; } + + // members + + std::vector minuit_internal_x_; + std::vector _grad; +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_LikelihoodGradientClad diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx index 95e08c158e4c5..2562575abaf26 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx @@ -11,6 +11,7 @@ */ #include "RooFit/TestStatistics/LikelihoodGradientWrapper.h" +#include "LikelihoodGradientClad.h" #include "RooMinimizer.h" // including derived classes for factory method @@ -71,6 +72,11 @@ LikelihoodGradientWrapper::create(LikelihoodGradientMode likelihoodGradientMode, RooMinimizer *minimizer) { switch (likelihoodGradientMode) { + case LikelihoodGradientMode::clad: { + return std::make_unique(std::move(likelihood), std::move(calculationIsClean), nDim, + minimizer); + break; + } case LikelihoodGradientMode::multiprocess: { #ifdef R__HAS_ROOFIT_MULTIPROCESS return std::make_unique(std::move(likelihood), std::move(calculationIsClean), nDim,