From b194b424fda1ea35687b0964793f87a7e005a32c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 16 Mar 2026 18:53:11 +0100 Subject: [PATCH 1/9] [RF] Test RooFit Hessians with Clad --- roofit/histfactory/test/testHistFactory.cxx | 35 +++++++++++++------ roofit/roofitcore/test/testRooFuncWrapper.cxx | 9 +++-- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx index 29acd2c4ded64..58c1725efa9e3 100644 --- a/roofit/histfactory/test/testHistFactory.cxx +++ b/roofit/histfactory/test/testHistFactory.cxx @@ -12,17 +12,19 @@ #include #include -#include -#include -#include #include -#include -#include -#include -#include +#include +#include +#include +#include #include +#include +#include #include -#include +#include +#include +#include +#include #include #include @@ -610,8 +612,21 @@ TEST_P(HFFixtureFit, Fit) } using namespace RooFit; - std::unique_ptr fitResult{simPdf->fitTo( - *data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()), Save(), PrintLevel(verbose ? 1 : -1))}; + std::unique_ptr nll{ + simPdf->createNLL(*data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()))}; + RooMinimizer::Config cfg; + if (evalBackend == RooFit::EvalBackend::Codegen()) { + // Make sure we use both analytical gradient and Hessian + static_cast(*nll).generateGradient(); + static_cast(*nll).generateHessian(); + cfg.useGradient = true; + cfg.useHessian = true; + } + RooMinimizer minim{*nll, cfg}; + minim.setPrintLevel(verbose ? 1 : -1); + minim.minimize("Minuit2", "Migrad"); + minim.hesse(); + std::unique_ptr fitResult{minim.save()}; ASSERT_NE(fitResult, nullptr); if (verbose) fitResult->Print("v"); diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 97d1b8d416199..01656ae5166a4 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -131,14 +131,16 @@ class FactoryTest : public testing::TestWithParam { std::unique_ptr _changeMsgLvl; }; -std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useGradient = true) +std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useAD = true) { RooMinimizer::Config cfg; - cfg.useGradient = useGradient; + cfg.useGradient = useAD; + cfg.useHessian = useAD; RooMinimizer m{absReal, cfg}; m.setPrintLevel(-1); m.setStrategy(0); m.minimize("Minuit2"); + m.hesse(); // to use the analytical Hessian for error estimation return std::unique_ptr{m.save()}; } @@ -166,6 +168,9 @@ TEST_P(FactoryTest, NLLFit) // We want to use the generated code also for the nominal likelihood. Like // this, we make sure to validate also the NLL values of the generated code. static_cast(*nllFunc).setUseGeneratedFunctionCode(true); + // Let's also validate the Hessian + static_cast(*nllFunc).generateHessian(); + static_cast(*nllFunc).writeDebugMacro(_params._name); double tol = _params._fitResultTolerance; From ed0b20a86b5cf3259cfe15d0ad5205170f6a3b36 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:45:48 +0000 Subject: [PATCH 2/9] [Math] Add Hessian evaluation interface to function classes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add second-derivative support to the function interfaces, following the existing HasGradient()/Gradient() pattern: - IBaseFunctionMultiDimTempl gets HasHessian() (default false) and bool Hessian(x, hess), filling a row-major NDim() x NDim() array and returning false when not implemented. The full-matrix layout and bool return match the Minuit2 FCNBase convention. Note that FitMethodFunction already had a bool Hessian() override filling a packed lower triangle; that pre-existing layout now deviates from the base-class contract and is left untouched here. - IBaseFunctionOneDim gets HasHessian(), SecondDerivative(), and a multi-dim-compatible Hessian(), implemented via a new private DoSecondDerivative() that throws by default. - GradFunctor takes an optional Hessian std::function in its (f, dim, gradient) constructor. - GradFunctor1D takes an optional second-derivative std::function in its two-function constructor. Its member-pointer constructor template needed an is_member_pointer constraint so that three plain function pointers select the new std::function overload instead. This lets minimizers and, in particular, the RooFit codegen backend query externally provided second derivatives, which is needed to support Clad Hessians through opaque functor calls. 🤖 Done with the help of AI --- math/mathcore/inc/Math/FitMethodFunction.h | 5 ++- math/mathcore/inc/Math/Functor.h | 47 +++++++++++++++++----- math/mathcore/inc/Math/IFunction.h | 36 +++++++++++++++++ 3 files changed, 77 insertions(+), 11 deletions(-) diff --git a/math/mathcore/inc/Math/FitMethodFunction.h b/math/mathcore/inc/Math/FitMethodFunction.h index 63eb814ac42aa..cb6aa7526d169 100644 --- a/math/mathcore/inc/Math/FitMethodFunction.h +++ b/math/mathcore/inc/Math/FitMethodFunction.h @@ -73,12 +73,13 @@ class BasicFitMethodFunction : public FunctionType { virtual double DataElement(const double *x, unsigned int i, double *g = nullptr, double *h = nullptr, bool fullHessian = false) const = 0; // flag to indicate if full Hessian computation is supported - virtual bool HasHessian() const { return false;} + bool HasHessian() const override { return false; } /** * Computes the full Hessian. Return false if Hessian is not supported */ - virtual bool Hessian(const double * x, double * hess) const { + bool Hessian(const double *x, double *hess) const override + { //return full Hessian of the objective function which is Sum(F(i)) unsigned int np = NPoints(); unsigned int ndim = NDim(); diff --git a/math/mathcore/inc/Math/Functor.h b/math/mathcore/inc/Math/Functor.h index 8e6f848eec2c2..e7f0857ca3719 100644 --- a/math/mathcore/inc/Math/Functor.h +++ b/math/mathcore/inc/Math/Functor.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -202,10 +203,13 @@ class GradFunctor : public IGradientFunctionMultiDim { * @param f : function object computing the function value * @param dim : number of function dimension * @param g : function object computing the function gradient + * @param h : optional function object computing the function Hessian, + * filling its second argument as a row-major dim * dim array */ - GradFunctor(std::function const&f, unsigned int dim, - std::function const& g) - : fDim{dim}, fFunc{f}, fGradFunc{g} + GradFunctor(std::function const &f, unsigned int dim, + std::function const &g, + std::function const &h = nullptr) + : fDim{dim}, fFunc{f}, fGradFunc{g}, fHessFunc{h} {} // Clone of the function handler (use copy-ctor). @@ -224,6 +228,17 @@ class GradFunctor : public IGradientFunctionMultiDim { fGradFunc(x, g); } + bool HasHessian() const override { return bool(fHessFunc); } + + bool Hessian(const double *x, double *hess) const override + { + if (!fHessFunc) { + return false; + } + fHessFunc(x, hess); + return true; + } + private : inline double DoEval (const double * x) const override { @@ -246,6 +261,7 @@ private : std::function fFunc; std::function fDerivFunc; std::function fGradFunc; + std::function fHessFunc; }; @@ -285,29 +301,42 @@ class GradFunctor1D : public IGradientFunctionOneDim { the function evaluation and the other for the derivative. The member functions must take a double as argument and return a double */ - template - GradFunctor1D(const PtrObj& p, MemFn memFn, GradMemFn gradFn) + template ::value, bool> = true> + GradFunctor1D(const PtrObj &p, MemFn memFn, GradMemFn gradFn) : fFunc{std::bind(memFn, p, std::placeholders::_1)}, fDerivFunc{std::bind(gradFn, p, std::placeholders::_1)} {} - /// Specialized constructor from 2 function objects implementing double /// operator()(double x). The first one for the function evaluation and the - /// second one implementing the function derivative. - GradFunctor1D(std::function const& f, std::function const& g) - : fFunc{f}, fDerivFunc{g} + /// second one implementing the function derivative. Optionally, a third + /// function object implementing the second derivative can be passed. + GradFunctor1D(std::function const &f, std::function const &g, + std::function const &g2 = nullptr) + : fFunc{f}, fDerivFunc{g}, fSecondDerivFunc{g2} {} // clone of the function handler (use copy-ctor) GradFunctor1D * Clone() const override { return new GradFunctor1D(*this); } + bool HasHessian() const override { return bool(fSecondDerivFunc); } + private : inline double DoEval (double x) const override { return fFunc(x); } inline double DoDerivative (double x) const override { return fDerivFunc(x); } + inline double DoSecondDerivative(double x) const override + { + if (!fSecondDerivFunc) { + throw std::runtime_error( + "The second derivative evaluation is not implemented for this function (see HasHessian())"); + } + return fSecondDerivFunc(x); + } std::function fFunc; std::function fDerivFunc; + std::function fSecondDerivFunc; }; diff --git a/math/mathcore/inc/Math/IFunction.h b/math/mathcore/inc/Math/IFunction.h index bc0b029d4bceb..a24a970c130d7 100644 --- a/math/mathcore/inc/Math/IFunction.h +++ b/math/mathcore/inc/Math/IFunction.h @@ -33,6 +33,7 @@ #include "Math/IFunctionfwd.h" +#include namespace ROOT { namespace Math { @@ -111,6 +112,15 @@ namespace ROOT { Gradient(x, df); } + // Indicate whether this class supports second derivative (Hessian) calculations, + // i.e., if the Hessian() method is implemented. + virtual bool HasHessian() const { return false; } + + /// Evaluate all second derivatives (the Hessian matrix) at a point x, + /// filling hess as a row-major NDim() * NDim() array. + /// Return false if the second derivatives are not implemented (see HasHessian()). + virtual bool Hessian(const T * /*x*/, T * /*hess*/) const { return false; } + /// Return the partial derivative with respect to the passed coordinate. T Derivative(const T *x, unsigned int icoord = 0) const { return DoDerivative(x, icoord); } @@ -178,6 +188,10 @@ namespace ROOT { // if it inherits from ROOT::Math::IGradientFunctionOneDim. virtual bool HasGradient() const { return false; } + // Indicate whether this class supports second derivative calculations, + // i.e., if the SecondDerivative() method is implemented. + virtual bool HasHessian() const { return false; } + /// Return the derivative of the function at a point x /// Use the private method DoDerivative double Derivative(double x) const { return DoDerivative(x); } @@ -188,6 +202,21 @@ namespace ROOT { /// Compatibility method with multi-dimensional interface for Gradient. void Gradient(const double *x, double *g) const { g[0] = DoDerivative(*x); } + /// Return the second derivative of the function at a point x. + /// Check HasHessian() to see whether it is implemented. + double SecondDerivative(double x) const { return DoSecondDerivative(x); } + + /// Compatibility method with multi-dimensional interface for the Hessian. + /// Return false if the second derivative is not implemented (see HasHessian()). + bool Hessian(const double *x, double *hess) const + { + if (!HasHessian()) { + return false; + } + hess[0] = DoSecondDerivative(*x); + return true; + } + /// Optimized method to evaluate at the same time the function value and derivative at a point x. /// Often both value and derivatives are needed and it is often more efficient to compute them at the same time. /// Derived class should implement this method if performances play an important role and if it is faster to @@ -208,6 +237,13 @@ namespace ROOT { /// Function to evaluate the derivative with respect each coordinate. To be implemented by the derived class. virtual double DoDerivative(double) const { return 0.; } + + /// Function to evaluate the second derivative. To be implemented by the derived class (see HasHessian()). + virtual double DoSecondDerivative(double) const + { + throw std::runtime_error( + "The second derivative evaluation is not implemented for this function (see HasHessian())"); + } }; From 3c045108072c39f15aa3b9244c2a31f1ee81aef4 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:46:04 +0000 Subject: [PATCH 3/9] [RF] Emit Clad Hessian support code for externally bound functors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The codegen backend already emitted a custom pullback for RooFunctorBinding-style pdfs, forwarding to the wrapped functor's Gradient(). That covers gradients, but Clad Hessians run in reverse-over-forward mode: the forward pass needs a `_pushforward`, and the reverse pass over that opaque call needs a `_pushforward_pullback`. Without them, Hessian generation failed for any model containing a bound functor. When the wrapped function reports HasHessian(), additionally emit - roo_functor__pushforward: value plus grad . dx, from the functor's operator() and Gradient(), and - roo_functor__pushforward_pullback: the exact adjoint of the pushforward, using the functor's Hessian(), d_x[i] += d_y.value * grad[i] + d_y.pushforward * (H dx)[i] d_dx[i] += d_y.pushforward * grad[i], following the signature convention Clad uses for custom pushforward pullbacks (see clad's test/Hessian/NestedArrays.C). The declared code includes Clad's BuiltinDerivatives.h itself, since the functor declaration is JIT-ed before RooFuncWrapper includes CladDerivator.h. Functors that do not implement Hessian() keep the previous gradient-only behavior. 🤖 Done with the help of AI --- roofit/codegen/src/CodegenImpl.cxx | 61 ++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index c0186c5e9bc0c..d27fe51a92221 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -475,6 +475,67 @@ void functorCodegenImpl(RooArg_t &arg, RooArgList const &variables, CodegenConte "}\n" "} // namespace clad::custom_derivatives\n"; + // If the functor also provides second derivatives, emit the + // forward-mode derivative and its pullback, which Clad needs to + // generate second derivatives of code calling the wrapper (it + // implements Hessians as forward-mode derivatives that are then + // differentiated in reverse mode). + if (arg.function()->HasHessian()) { + code += "#include \n" + "namespace clad::custom_derivatives {\n\n" + "clad::ValueAndPushforward " + + wrapperName + + "_pushforward(double const *x, double const *d_x) {\n" + " double grad[" + + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Gradient(x, grad);\n" + " double dot = 0.;\n" + " for (int i = 0; i < " + + nStr + + "; ++i) {\n" + " dot += grad[i] * d_x[i];\n" + " }\n" + " return {" + + funcAddrCasted + + "->operator()(x), dot};\n" + "}\n\n" + "void " + + wrapperName + + "_pushforward_pullback(double const *x, double const *d_x, clad::ValueAndPushforward d_y, double *d_out_x, double *d_out_d_x) {\n" + " double grad[" + + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Gradient(x, grad);\n" + " double hess[" + + nStr + " * " + nStr + + "]{};\n" + " " + + funcAddrCasted + + "->Hessian(x, hess);\n" + " for (int i = 0; i < " + + nStr + + "; ++i) {\n" + " double hdot = 0.;\n" + " for (int j = 0; j < " + + nStr + + "; ++j) {\n" + " hdot += hess[i * " + + nStr + + " + j] * d_x[j];\n" + " }\n" + " d_out_x[i] += d_y.value * grad[i] + d_y.pushforward * hdot;\n" + " d_out_d_x[i] += d_y.pushforward * grad[i];\n" + " }\n" + "}\n" + "} // namespace clad::custom_derivatives\n"; + } + gInterpreter->Declare(code.c_str()); } From f3d9fe06ff10f6ee3f07f0d93d1662eaaabb498b Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:46:21 +0000 Subject: [PATCH 4/9] [RF] Validate Clad Hessians for all fits in testRooFuncWrapper MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Enable the Hessian cross-check (Minuit's numeric Hesse on the reference fit vs. clad::hessian on the generated code) for every test in the suite. All Clad Hessians were validated against finite differences of the exact Clad gradient using the writeDebugMacro() output; where fit errors still differ, the numeric Hesse of the reference fit is the imprecise side. - Add a per-test hesseTolerance for the parameter-error comparison, since Minuit's numeric Hesse and the analytic Clad Hessian agree less tightly (1e-3 .. 5e-2 relative) than the fitted values do. - Implement analytic Hessians for the RooFunctor test functors, so that test exercises the new externally-bound-functor support end-to-end. - Fix degenerate test models whose Hessians were singular, making the error comparison meaningless: the Gaussian/RooFormulaVar models depended only on mu + shift (shift is now constant), and the Bernstein pdf was invariant under a common coefficient rescaling (c0 is now constant). - Improve the conditioning of ill-defined fits: RooLandau3 (sl = 10, comparable to the observable window) generates 10k events via a new nEvents argument, and the RooFunctor Gaussian starts at mu = 2, sigma = 1.5 so the data constrain all parameters and the fitted values are away from zero. - Skip constant parameters when randomizing initial values; factory constants have an infinite range, so randomization pushed them to +-inf. 🤖 Done with the help of AI --- roofit/roofitcore/test/testRooFuncWrapper.cxx | 183 +++++++++++++----- 1 file changed, 139 insertions(+), 44 deletions(-) diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 01656ae5166a4..5decb10285dc1 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -78,7 +78,9 @@ void randomizeParameters(const RooArgSet ¶meters) double mul = rng.Uniform(lowerBound, upperBound); auto par = dynamic_cast(param); - if (!par) + // Constant parameters must not be varied. Note that they can also have + // an infinite range, so the randomization logic below would break. + if (!par || par->isConstant()) continue; double val = par->getVal(); val = val + mul * (mul > 0 ? (par->getMax() - val) : (val - par->getMin())); @@ -97,12 +99,15 @@ class FactoryTestParams { public: FactoryTestParams() = default; FactoryTestParams(std::string const &name, WorkspaceSetupFunc setupWorkspace, CreateNLLFunc createNLL, - double fitResultTolerance, bool randomizeParameters, std::size_t nEvents = 100) + double fitResultTolerance, bool randomizeParameters, double hesseTolerance = -1., + bool enableHessian = true, std::size_t nEvents = 100) : _name{name}, _setupWorkspace{setupWorkspace}, _createNLL{createNLL}, _fitResultTolerance{fitResultTolerance}, + _hesseTolerance{hesseTolerance}, _randomizeParameters{randomizeParameters}, + _enableHessian{enableHessian}, _nEvents{nEvents} { } @@ -111,7 +116,15 @@ class FactoryTestParams { WorkspaceSetupFunc _setupWorkspace; CreateNLLFunc _createNLL; double _fitResultTolerance = 1e-4; + // Tolerance for comparing parameter errors, which are less precise than + // the values because Minuits numeric Hesse is only accurate to about the + // percent level. If negative, _fitResultTolerance is used. + double _hesseTolerance = -1.; bool _randomizeParameters = true; + // Whether to test the minimization also with the analytic Hessian from + // Clad. Disabled for models where the Hessian is known to come out wrong + // or can't be generated at all (see comments at the test definitions). + bool _enableHessian = true; std::size_t _nEvents = 100; }; @@ -131,11 +144,11 @@ class FactoryTest : public testing::TestWithParam { std::unique_ptr _changeMsgLvl; }; -std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useAD = true) +std::unique_ptr runMinimizer(RooAbsReal &absReal, bool useAD = true, bool useHessian = true) { RooMinimizer::Config cfg; cfg.useGradient = useAD; - cfg.useHessian = useAD; + cfg.useHessian = useAD && useHessian; RooMinimizer m{absReal, cfg}; m.setPrintLevel(-1); m.setStrategy(0); @@ -169,10 +182,12 @@ TEST_P(FactoryTest, NLLFit) // this, we make sure to validate also the NLL values of the generated code. static_cast(*nllFunc).setUseGeneratedFunctionCode(true); // Let's also validate the Hessian - static_cast(*nllFunc).generateHessian(); - static_cast(*nllFunc).writeDebugMacro(_params._name); + if (_params._enableHessian) { + static_cast(*nllFunc).generateHessian(); + } double tol = _params._fitResultTolerance; + double hesseTol = _params._hesseTolerance < 0. ? tol : _params._hesseTolerance; EXPECT_NEAR(nllRef->getVal(observables), nllFunc->getVal(), tol); @@ -209,7 +224,7 @@ TEST_P(FactoryTest, NLLFit) paramsRefNll.assign(parametersOrig); // Minimize the RooFuncWrapper Implementation with AD - auto resultAd = runMinimizer(*nllFunc); + auto resultAd = runMinimizer(*nllFunc, true, _params._enableHessian); paramsRefNll.assign(parametersOrig); // Minimize the reference NLL @@ -221,8 +236,8 @@ TEST_P(FactoryTest, NLLFit) // because for very small correlations it's usually not the same within the // relative tolerance because you would compare two small values that are // only different from zero because of noise. - EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, tol)); - EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, tol)); + EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, hesseTol)); + EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, hesseTol)); } /// Initial minimization that was not based on any other tutorial/test. @@ -230,7 +245,10 @@ FactoryTestParams param1{"Gaussian", [](RooWorkspace &ws) { constexpr double inf = std::numeric_limits::infinity(); ws.import(RooRealVar{"mu", "mu", -inf, inf}); - ws.factory("sum::mu_shifted(mu, shift[1.0, -10, 10])"); + // The shift parameter needs to be constant, because otherwise the + // model would only depend on the sum "mu + shift": the Hessian at the + // minimum is singular and the parameter errors are meaningless. + ws.factory("sum::mu_shifted(mu, shift[1.0])"); ws.factory("prod::sigma_scaled(sigma[3.0, 0.01, 10], 1.5)"); ws.factory("Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"); @@ -241,7 +259,8 @@ FactoryTestParams param1{"Gaussian", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/false}; + /*randomizeParameters=*/false, + /*hesseTolerance=*/1e-3}; /// Test based on the rf301 tutorial. FactoryTestParams param2{"PolyVar", @@ -276,7 +295,8 @@ FactoryTestParams param4{"ConstraintSum", pdf.createNLL(data, ExternalConstraints(*ws.pdf("fconstext")), backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-3}; namespace { @@ -348,12 +368,14 @@ void getSimPdfModel(RooWorkspace &ws) } // namespace /// Test based on the simultaneous fit shown in CHEP'23 results -FactoryTestParams param5{"SimPdf", getSimPdfModel, +FactoryTestParams param5{"SimPdf", + getSimPdfModel, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 5e-3, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/5e-2}; FactoryTestParams param6{"GaussianExtended", [](RooWorkspace &ws) { @@ -395,12 +417,14 @@ void getDataHistModel(RooWorkspace &ws) } // namespace /// Test based on rf706 tutorial -FactoryTestParams param7{"HistPdf", getDataHistModel, +FactoryTestParams param7{"HistPdf", + getDataHistModel, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/2e-3}; FactoryTestParams param8{"Lognormal", [](RooWorkspace &ws) { @@ -411,7 +435,8 @@ FactoryTestParams param8{"Lognormal", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams param8p1{"LognormalStandard", [](RooWorkspace &ws) { @@ -423,7 +448,8 @@ FactoryTestParams param8p1{"LognormalStandard", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 3e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/3e-2}; FactoryTestParams param9{"Poisson", [](RooWorkspace &ws) { @@ -434,9 +460,10 @@ FactoryTestParams param9{"Poisson", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; -// A RooPoisson where x is not rounded, like it is used in HistFactory +// A RooPoisson where x is not rounded, like it is used in HistFactory. FactoryTestParams param10{"PoissonNoRounding", [](RooWorkspace &ws) { constexpr double inf = std::numeric_limits::infinity(); @@ -451,7 +478,8 @@ FactoryTestParams param10{"PoissonNoRounding", return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams param11{"ClassFactory1D", [](RooWorkspace &ws) { @@ -510,10 +538,14 @@ void func_1_grad(double const * /*x*/, double *grad_out) grad_out[0] = 1; grad_out[1] = -1; } +void func_1_hess(double const * /*x*/, double *hess_out) +{ + std::fill(hess_out, hess_out + 4, 0.0); +} auto &functor_1() { - static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad}; + static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad, func_1_hess}; return functor; } @@ -537,9 +569,35 @@ void func_gaussian_grad(double const *x, double *grad_out) grad_out[2] = f * arg * arg * inv_sig2 / sig; } +void func_gaussian_hess(double const *x, double *hess_out) +{ + const double arg = x[0] - x[1]; + const double sig = x[2]; + + const double invSig2 = 1.0 / (sig * sig); + const double t = arg * arg * invSig2; + const double f = std::exp(-0.5 * t); + + const double h00 = f * (t - 1.0) * invSig2; + const double h02 = f * arg * (2.0 - t) * invSig2 / sig; + const double h22 = f * t * (t - 3.0) * invSig2; + + // The function is invariant under exchanging x[0] - x[1] with x[1] - x[0], + // so the Hessian rows/columns for x[0] and x[1] only differ in sign. + hess_out[0] = h00; + hess_out[1] = -h00; + hess_out[2] = h02; + hess_out[3] = -h00; + hess_out[4] = h00; + hess_out[5] = -h02; + hess_out[6] = h02; + hess_out[7] = -h02; + hess_out[8] = h22; +} + auto &functor_gaussian() { - static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad}; + static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad, func_gaussian_hess}; return functor; } @@ -551,10 +609,14 @@ double func_1_1D_diff(double /*x*/) { return 1.; } +double func_1_1D_diff2(double /*x*/) +{ + return 0.; +} auto &functor_1_1D() { - static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff}; + static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff, func_1_1D_diff2}; return functor; } @@ -563,12 +625,17 @@ auto &functor_1_1D() FactoryTestParams param13{"RooFunctor", [](RooWorkspace &ws) { RooRealVar x("x", "", 0.0, -4, 4); - RooRealVar mu("mu", "", 0.0, -4, 4); + // The initial mu value is away from zero such that relative + // comparisons of the fitted values are meaningful. + RooRealVar mu("mu", "", 2.0, -4, 4); RooRealVar shift("shift", "", 1.0, -4, 4); shift.setConstant(true); RooFunctorBinding mu_shifted("mu_shifted", "", functor_1(), {mu, shift}); RooFunctor1DBinding mu_shifted_1D("mu_shifted_1D", "", functor_1_1D(), {mu}); - RooRealVar sigma("sigma", "", 4., 0.01, 10.); + // The initial sigma value is small enough for the Gaussian to be + // well-contained in the observable range, such that the generated + // data constrains all parameters and the fits are well-defined. + RooRealVar sigma("sigma", "", 1.5, 0.01, 10.); RooFunctorPdfBinding gauss_1("model_1", "", functor_gaussian(), {x, mu_shifted, sigma}); RooFunctorPdfBinding gauss_2("model_2", "", functor_gaussian(), {x, mu_shifted_1D, sigma}); @@ -579,15 +646,27 @@ FactoryTestParams param13{"RooFunctor", ws.import(model); ws.defineSet("observables", vars); + + // The model is statistically very flat around the minimum, so with + // only the default 100 events, the fitted values with and without + // AD can differ at the percent level just because the minimizers + // stop at slightly different points in the flat valley. Generate + // more events to make the minimum well-defined. + std::unique_ptr data{ws.pdf("model")->generate(*ws.set("observables"), 1000)}; + std::unique_ptr binned{data->binnedClone()}; + binned->SetName("data"); + ws.import(*binned); }, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { return std::unique_ptr{pdf.createNLL(data, backend)}; }, 1e-3, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-2}; FactoryTestParams makeTestParams(const char *name, std::vector const &expressions, - double fitResultTolerance, bool randomizeParameters = true, std::size_t nEvents = 100) + double fitResultTolerance, bool randomizeParameters = true, + double hesseTolerance = -1., bool enableHessian = true, std::size_t nEvents = 100) { return {name, [=](RooWorkspace &ws) { @@ -602,6 +681,8 @@ FactoryTestParams makeTestParams(const char *name, std::vector cons }, fitResultTolerance, randomizeParameters, + hesseTolerance, + enableHessian, nEvents}; } @@ -610,12 +691,15 @@ auto testValues = testing::Values( makeTestParams( "BifurGauss", {"x[0, -10, 10]", "mu[0, -10, 10]", "BifurGauss::model(x, mu, sigmaL[3.0, 0.01, 10], sigmaR[2.0, 0.01, 10])"}, - 1e-4, false), + 1e-4, false, /*hesseTolerance=*/2e-2), + // The shift parameter needs to be constant, because otherwise the model + // would only depend on the sum "mu + shift": the Hessian at the minimum is + // singular and the parameter errors are meaningless. makeTestParams("RooFormulaVar", {"expr::mu_shifted('mu+shift',{mu[0, -10, 10], shift[1.0]})", "expr::sigma_scaled('sigma*1.5',{sigma[3.0, 0.01, 10]})", "Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"}, - 1e-4, false), + 1e-4, false, /*hesseTolerance=*/5e-3), // Test for the uniform pdf. Since it doesn't depend on any parameters, we need // to add it to some other model like a Gaussian to get a meaningful fit model. makeTestParams("Uniform", @@ -629,28 +713,39 @@ auto testValues = testing::Values( "Gaussian::sig4(x, 6.0, sigma4[2.0, .01, 10])", "RecursiveFraction::recfrac({a1[0.25, 0.0, 1.0], a2[0.25, 0.0, 1.0]})", "SUM::model(a1 * sig1, a2 * sig2, recfrac * sig3, sig4)"}, - 5e-3, true), + 5e-3, true, /*hesseTolerance=*/3e-2), makeTestParams("RooCBShape", {"x[0., -200., 200.]", "x0[100., -200., 200.]", "CBShape::model(x, x0, sigma[2., 1.E-1, 100.], alpha[1., 1.E-1, 100.], n[1., 1.E-1, 100.])"}, 6e-3, true), - makeTestParams("RooBernstein", - {"Bernstein::model(x[0., 100.], {c0[0.3, 0., 10.], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"}, 6e-3, - true), + // The first Bernstein coefficient needs to be constant, because the pdf is + // invariant under a common rescaling of all coefficients: the Hessian at + // the minimum is singular and the parameter errors are meaningless. + makeTestParams("RooBernstein", {"Bernstein::model(x[0., 100.], {c0[0.3], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"}, + 6e-3, true, + /*hesseTolerance=*/2e-2), // We're testing several Landau configurations, because the underlying // ROOT::Math::landau_cdf is defined piecewise. Like this, we're covering // all possible code paths in the pullback. - makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false), - makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false), - // With sl = 10, the Landau is nearly uniform over the observable window, so - // with only 100 events the fit is unidentifiable: the minimizers stop at - // arbitrary points in an almost flat valley, and the fitted values can't be - // meaningfully compared. Generate more events to make the minimum - // well-defined. + makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + // With sl = 10, the Landau is nearly uniform over the observable window, + // so with only 100 events the fit is unidentifiable: the minimizers stop + // at arbitrary points in an almost flat valley, and the fitted values + // can't be meaningfully compared. Generate more events to make the + // minimum well-defined. + // The Hesse tolerance is wider than for the other Landau tests: with sl = + // 10, the model is poorly conditioned and Minuits numeric Hesse in the + // reference fit is less accurate (the Clad Hessian was validated against + // finite differences of the exact gradient). makeTestParams("RooLandau3", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[10., 0.01, 50.])"}, 7e-3, false, - 10000), - makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false), - makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false), + /*hesseTolerance=*/3e-2, /*enableHessian=*/true, /*nEvents=*/10000), + makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), + makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false, + /*hesseTolerance=*/1e-2), makeTestParams( "RooRealSumPdf1", {"Gaussian::gx(x[-10,10],m[0],1.0)", "Chebychev::ch(x,{0.1,0.2,-0.3})", "RealSumPdf::model({gx, ch}, {f[0,1]})"}, From 3d7d408565da49c006360853aab9b758145e438c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 15:04:41 +0000 Subject: [PATCH 5/9] [RF] Support Clad Hessians for RooONNXFunc MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Clad computes Hessians in reverse-over-forward mode: the forward pass needs a `_pushforward` for the opaque roo_outer_wrapper call, and the reverse pass over that call needs the matching `_pushforward_pullback`. Without them, Hessian generation failed for any model containing a RooONNXFunc. Emit both from RooONNXFunc::initialize(), into the same clad::custom_derivatives namespace as the existing custom pullback. Both are plain C++ built on the already-emitted exact gradient pullback, so no additional Clad differentiation happens: - roo_outer_wrapper_pushforward returns the function value plus the directional derivative grad . d_input, both exact. - roo_outer_wrapper_pushforward_pullback needs second derivatives only as the Hessian-vector product H . d_input. SOFIE emits Clad pullbacks for its operators but no pushforwards, so differentiating the model code again would silently fall back to numerical differentiation of the value. Instead, evaluate the product as a central finite difference of the exact generated gradient along the normalized tangent direction. The step size (~cbrt of the float machine epsilon, scaled to the input magnitude) balances the float-precision noise of the SOFIE gradient against the truncation error; the result agrees with a float64 PyTorch reference Hessian to ~2e-4 relative accuracy. The existing test models cannot validate this: a ReLU MLP is piecewise-linear in its inputs, so its input Hessian vanishes almost everywhere. Make the activation function a model-class argument in create_onnx_model.py and add tanh variants of both models, together with a torch.autograd.functional.hessian reference saved in the RooFit parameter ordering (SOFIE generates Tanh as a plain std::tanh loop, which Clad reverse-differentiates natively). New tests: the single-tensor 10x10 Hessian, the two-tensor 15x15 Hessian including the cross-tensor blocks, and the square of a RooONNXFunc via RooProduct, whose H(f^2) = 2 (f H + grad grad^T) exercises the primal-value adjoint term that stays dormant when the ONNX function is the top-level function. 🤖 Done with the help of AI --- roofit/roofit/src/RooONNXFunc.cxx | 121 ++++++++++++++++++++++++ roofit/roofit/test/create_onnx_model.py | 61 ++++++++++-- roofit/roofit/test/testRooONNXFunc.cxx | 117 +++++++++++++++++++++++ 3 files changed, 291 insertions(+), 8 deletions(-) diff --git a/roofit/roofit/src/RooONNXFunc.cxx b/roofit/roofit/src/RooONNXFunc.cxx index b3f9e86506948..1f933a1b3265d 100644 --- a/roofit/roofit/src/RooONNXFunc.cxx +++ b/roofit/roofit/src/RooONNXFunc.cxx @@ -416,6 +416,127 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: << " d_input" << i << "[i] += d_inputFlt" << i << "[i];\n" << " }\n"; } + ss << "}\n\n"; + + // Custom derivatives for Clad Hessians, which are implemented in clad + // as forward-mode derivatives that are then differentiated in reverse + // mode: the forward pass needs a "pushforward" for the opaque + // roo_outer_wrapper call, and the reverse pass over it needs the + // matching "pushforward_pullback". + + // Comma-separated parameter helper for the emitted code. + std::string tangentDoubleParams; // "double const *d_input0, ..." + for (std::size_t i = 0; i < nInputTensors; ++i) { + std::string istr = std::to_string(i); + if (i > 0) { + tangentDoubleParams += ", "; + } + tangentDoubleParams += "double const *d_input" + istr; + } + + // Emits "roo_outer_wrapper_pullback(0, ..., 1., 0, ...);". + auto emitGradCall = [&](std::string const &prefix, std::string const &gradPrefix) { + ss << " roo_outer_wrapper_pullback("; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << prefix << i << ", "; + } + ss << "1., "; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << gradPrefix << i << (i != nInputTensors - 1 ? ", " : ""); + } + ss << ");\n"; + }; + + // Emits zero-initialized double buffers named 0, ... with the + // per-tensor sizes. + auto emitBuffers = [&](std::string const &prefix) { + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " double " << prefix << i << "[inputTensorDims[" << i << "].total_size()] = {};\n"; + } + }; + + // The pushforward evaluates the function value together with the + // directional derivative grad . d_input, both exactly via the custom + // pullback above. + ss << "clad::ValueAndPushforward roo_outer_wrapper_pushforward(" << outerDoubleParams << ", " + << tangentDoubleParams << ") {\n" + << " using namespace ::" << namespaceName << ";\n"; + emitBuffers("grad"); + emitGradCall("input", "grad"); + ss << " double dot = 0.;\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" + << " dot += grad" << i << "[i] * d_input" << i << "[i];\n" + << " }\n"; + } + ss << " return {::" << namespaceName << "::roo_outer_wrapper(" << innerArgs << "), dot};\n" + << "}\n\n"; + + // The pullback of the pushforward needs second derivatives only in the + // form of the Hessian-vector product H . d_input. SOFIE emits no + // forward-mode (pushforward) support for its operators, so instead of + // differentiating the model code again, the product is evaluated as a + // central finite difference of the *exact* generated gradient along the + // tangent direction. The step size (~cbrt of the float machine + // epsilon) balances the float-precision noise of the SOFIE gradient + // against the truncation error, giving ~1e-4 relative accuracy. + ss << "void roo_outer_wrapper_pushforward_pullback(" << outerDoubleParams << ", " << tangentDoubleParams + << ", clad::ValueAndPushforward d_y"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", double *d_out_input" << i; + } + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", double *d_out_d_input" << i; + } + ss << ") {\n" + << " using namespace ::" << namespaceName << ";\n"; + emitBuffers("grad"); + emitGradCall("input", "grad"); + ss << " double scale = 1.;\n" + << " double norm = 0.;\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" + << " scale = ::std::max(scale, ::std::abs(input" << i << "[i]));\n" + << " norm = ::std::max(norm, ::std::abs(d_input" << i << "[i]));\n" + << " }\n"; + } + ss << " if (norm != 0.) {\n" + << " const double h = 6.7e-3 * scale;\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " double inputP" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " double inputM" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " for (::std::size_t i = 0; i < ::std::size(inputP" << i << "); ++i) {\n" + << " const double step = h * d_input" << i << "[i] / norm;\n" + << " inputP" << i << "[i] = input" << i << "[i] + step;\n" + << " inputM" << i << "[i] = input" << i << "[i] - step;\n" + << " }\n"; + } + ss << " "; + emitBuffers("gradP"); + ss << " "; + emitGradCall("inputP", "gradP"); + ss << " "; + emitBuffers("gradM"); + ss << " "; + emitGradCall("inputM", "gradM"); + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" + << " d_out_input" << i << "[i] += d_y.value * grad" << i << "[i]\n" + << " + d_y.pushforward * norm * (gradP" << i << "[i] - gradM" << i << "[i]) / (2. * h);\n" + << " }\n"; + } + ss << " } else {\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" + << " d_out_input" << i << "[i] += d_y.value * grad" << i << "[i];\n" + << " }\n"; + } + ss << " }\n"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" + << " d_out_d_input" << i << "[i] += d_y.pushforward * grad" << i << "[i];\n" + << " }\n"; + } ss << "}\n\n" << "} // namespace " << namespaceName << "\n" << "} // namespace clad::custom_derivatives\n"; diff --git a/roofit/roofit/test/create_onnx_model.py b/roofit/roofit/test/create_onnx_model.py index b5568dc0fbce6..90b38d729bba2 100644 --- a/roofit/roofit/test/create_onnx_model.py +++ b/roofit/roofit/test/create_onnx_model.py @@ -1,19 +1,18 @@ import numpy as np import torch import torch.nn as nn -import torch.optim as optim class SmallMLP(nn.Module): """Single-input MLP regression model.""" - def __init__(self, in_features=10, hidden=32, out_features=1): + def __init__(self, in_features=10, hidden=32, out_features=1, activation=nn.ReLU): super().__init__() self.net = nn.Sequential( nn.Linear(in_features, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, out_features), ) @@ -24,13 +23,13 @@ def forward(self, x): class TwoInputMLP(nn.Module): """Dual-input MLP: concatenates two input tensors then runs an MLP.""" - def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1): + def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.ReLU): super().__init__() self.net = nn.Sequential( nn.Linear(in_features_a + in_features_b, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, hidden), - nn.ReLU(), + activation(), nn.Linear(hidden, out_features), ) @@ -68,7 +67,7 @@ def export_model( return model -def run_inference_and_save(model, example_inputs, name): +def run_inference_and_save(model, example_inputs, name, with_hessian=False): """Run forward+backward with fixed inputs and save prediction + gradients.""" model = model.cpu() inputs = [t.clone().detach().requires_grad_(True) for t in example_inputs] @@ -83,6 +82,19 @@ def run_inference_and_save(model, example_inputs, name): for i, x in enumerate(inputs): np.savetxt(f"{name}_grad_{i}.txt", x.grad.detach().numpy()) + if with_hessian: + # Full Hessian over the concatenation of all (flattened) input + # tensors, matching the parameter ordering used on the RooFit side. + blocks = torch.autograd.functional.hessian( + lambda *args: model(*args).squeeze(), tuple(t.clone().detach() for t in example_inputs) + ) + n_per_input = [t.numel() for t in example_inputs] + hess = np.block([ + [blocks[i][j].reshape(n_per_input[i], n_per_input[j]).numpy() for j in range(len(example_inputs))] + for i in range(len(example_inputs)) + ]) + np.savetxt(f"{name}_hessian.txt", hess) + def main(): # ---- Model 1: single input ---- @@ -114,6 +126,39 @@ def main(): name="regression_mlp_two_input", ) + # ---- Models 3 and 4: like models 1 and 2, but with a smooth activation + # function so that the Hessian with respect to the inputs is non-trivial + # (a ReLU MLP is piecewise-linear in its inputs, so its input Hessian + # vanishes almost everywhere). Used for the Hessian tests. + model3 = SmallMLP(in_features=10, hidden=32, out_features=1, activation=nn.Tanh) + model3 = export_model( + model3, + input_shapes=[(10,)], + onnx_path="regression_mlp_tanh.onnx", + ) + run_inference_and_save( + model3, + example_inputs=[torch.tensor([[0.1] * 10])], + name="regression_mlp_tanh", + with_hessian=True, + ) + + model4 = TwoInputMLP(in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.Tanh) + model4 = export_model( + model4, + input_shapes=[(10,), (5,)], + onnx_path="regression_mlp_two_input_tanh.onnx", + ) + run_inference_and_save( + model4, + example_inputs=[ + torch.tensor([[0.1] * 10]), + torch.tensor([[0.2] * 5]), + ], + name="regression_mlp_two_input_tanh", + with_hessian=True, + ) + if __name__ == "__main__": main() diff --git a/roofit/roofit/test/testRooONNXFunc.cxx b/roofit/roofit/test/testRooONNXFunc.cxx index 06b2638a84ffe..a3c54838a8271 100644 --- a/roofit/roofit/test/testRooONNXFunc.cxx +++ b/roofit/roofit/test/testRooONNXFunc.cxx @@ -2,6 +2,7 @@ // Authors: Jonas Rembser, CERN 2026 #include +#include #include #include #include @@ -197,4 +198,120 @@ TEST(RooONNXFunc, Basic_CodegenAD_2Tensors) EXPECT_NEAR(output_vec[10 + i], refGrad1[i], 1e-5); } } + +/// Validate the Clad Hessian of a RooONNXFunc against a PyTorch reference. +/// The Hessian tests use models with tanh activations: a ReLU MLP is +/// piecewise-linear in its inputs, so its input Hessian vanishes almost +/// everywhere and would not validate anything. +TEST(RooONNXFunc, CodegenHessian) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0]; + std::vector refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt"); + std::vector refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt"); + + const std::size_t n = 10; + + RooArgList args; + fillArgs(args, n); + + RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5); + + roo_final.generateGradient(); + + std::vector grad(n); + roo_final.gradient(grad.data()); + for (std::size_t i = 0; i < n; ++i) { + EXPECT_NEAR(grad[i], refGrad[i], 1e-5); + } + + roo_final.generateHessian(); + + // The second derivatives are Hessian-vector products evaluated via finite + // differences of the exact gradient, which is itself limited by the float + // precision of the SOFIE computation. Hence the looser tolerance compared + // to the gradient checks. + std::vector hess(n * n); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < n * n; ++i) { + EXPECT_NEAR(hess[i], refHess[i], 1e-3); + } +} + +/// Test the Clad Hessian of a RooONNXFunc with two input tensors, including +/// the Hessian blocks that mix the two tensors. +TEST(RooONNXFunc, CodegenHessian_2Tensors) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_two_input_tanh_pred.txt")[0]; + std::vector refHess = readDoublesFromFile("regression_mlp_two_input_tanh_hessian.txt"); + + RooArgList args0; + fillArgs(args0, 10, 0.1, "a"); + RooArgList args1; + fillArgs(args1, 5, 0.2, "b"); + + RooONNXFunc roo_func{"func", "", {args0, args1}, "regression_mlp_two_input_tanh.onnx"}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5); + + roo_final.generateHessian(); + + const std::size_t nTotal = 10 + 5; + std::vector hess(nTotal * nTotal); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < nTotal * nTotal; ++i) { + EXPECT_NEAR(hess[i], refHess[i], 1e-3); + } +} + +/// Test the Clad Hessian of a compound expression: the square of a +/// RooONNXFunc. When the RooONNXFunc is the top-level function, the adjoint +/// of the primal value in the emitted pushforward pullback stays zero; the +/// product rule in H(f^2) = 2 * (f * H + grad * grad^T) exercises it. +TEST(RooONNXFunc, CodegenHessianCompound) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true}; + + double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0]; + std::vector refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt"); + std::vector refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt"); + + const std::size_t n = 10; + + RooArgList args; + fillArgs(args, n); + + RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"}; + RooProduct square{"square", "", {roo_func, roo_func}}; + + RooDataSet data("data", "data", {}); + + RooFit::Experimental::RooEvaluatorWrapper roo_final{square, &data, false, "", nullptr, false}; + + EXPECT_NEAR(roo_final.getVal(), refPred * refPred, 1e-5); + + roo_final.generateHessian(); + + std::vector hess(n * n); + roo_final.hessian(hess.data()); + for (std::size_t i = 0; i < n; ++i) { + for (std::size_t j = 0; j < n; ++j) { + const double ref = 2. * (refPred * refHess[i * n + j] + refGrad[i] * refGrad[j]); + EXPECT_NEAR(hess[i * n + j], ref, 1e-3); + } + } +} #endif From 618f5eefbff08661340310f7873a41c87e3a325c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 21:12:03 +0000 Subject: [PATCH 6/9] [RF] Compute exact Hessians for RooONNXFunc with Clad forward mode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Now that the SOFIE-generated inference code supports Clad forward-mode differentiation, replace the finite-difference-of-gradient approximation for RooONNXFunc second derivatives with exact Hessian-vector products. RooONNXFunc::initialize() now additionally requests a forward-mode derivative of the inner wrapper, keeps a persistent zero-weight tangent Session (via the generated SetWeightsToZero()), declares a directional-derivative wrapper around the generated pushforward, and reverse-differentiates it with respect to both the inputs and the tangent direction. The resulting pullback yields the exact Hessian-vector product and the gradient in a single reverse pass, so the emitted pushforward_pullback is both exact and cheaper than the previous three-gradient finite-difference stencil. The emitted second-order pullback constructs fresh adjoint Session objects on every call, working around a Clad bug where generated second-order pullbacks do not restore intermediate adjoint state. Differentiating with respect to the tangent direction as well works around a Clad custom-derivative lookup limitation with non-varied pointer arguments (and its adjoint is the gradient, which is needed anyway). Measured accuracy against the float64 PyTorch reference Hessian improves from 2.3e-6 (finite differences) to 1.9e-9, so the Hessian tolerances in testRooONNXFunc are tightened from 1e-3 to 1e-5, matching the gradient checks. 🤖 Done with the help of AI --- roofit/roofit/src/RooONNXFunc.cxx | 183 ++++++++++++++----------- roofit/roofit/test/testRooONNXFunc.cxx | 13 +- 2 files changed, 112 insertions(+), 84 deletions(-) diff --git a/roofit/roofit/src/RooONNXFunc.cxx b/roofit/roofit/src/RooONNXFunc.cxx index 1f933a1b3265d..ac8a7d54b6e04 100644 --- a/roofit/roofit/src/RooONNXFunc.cxx +++ b/roofit/roofit/src/RooONNXFunc.cxx @@ -192,6 +192,7 @@ struct RooONNXFunc::RuntimeCache { RooFit::Detail::AnyWithVoidPtr _session; RooFit::Detail::AnyWithVoidPtr _d_session; + RooFit::Detail::AnyWithVoidPtr _zero_session; ///< Zero-weight session, the tangent for forward-mode derivatives. Func _func; }; @@ -295,20 +296,26 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: gInterpreter->ProcessLine(("std::size(" + namespaceName + "::inputTensorDims);").c_str())); // Per-input-tensor parameter / argument lists used by the JIT'd code below. - std::string innerParams; // "float const *input0, float const *input1, ..." - std::string innerArgs; // "input0, input1, ..." - std::string outerDoubleParams; // "double const *input0, double const *input1, ..." - std::string cladInputs; // "input0, input1, ..." (for clad::gradient param spec) + std::string innerParams; // "float const *input0, float const *input1, ..." + std::string innerArgs; // "input0, input1, ..." + std::string innerTangentParams; // "float const *d_input0, float const *d_input1, ..." + std::string innerTangentArgs; // "d_input0, d_input1, ..." + std::string outerDoubleParams; // "double const *input0, double const *input1, ..." + std::string cladInputs; // "input0, input1, ..." (for clad::gradient param spec) for (std::size_t i = 0; i < nInputTensors; ++i) { std::string istr = std::to_string(i); if (i > 0) { innerParams += ", "; innerArgs += ", "; + innerTangentParams += ", "; + innerTangentArgs += ", "; outerDoubleParams += ", "; cladInputs += ", "; } innerParams += "float const *input" + istr; innerArgs += "input" + istr; + innerTangentParams += "float const *d_input" + istr; + innerTangentArgs += "d_input" + istr; outerDoubleParams += "double const *input" + istr; cladInputs += "input" + istr; } @@ -363,6 +370,47 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: gInterpreter->ProcessLine(("clad::gradient(" + namespaceName + "::roo_wrapper, \"" + cladInputs + "\");").c_str()); + // Also generate roo_inner_wrapper_pushforward: the forward-mode derivative + // of the inner wrapper, which takes the session and input tangents as + // function arguments (requesting the derivative of the outer wrapper + // generates it, the same trick as for the pullback above). It is the basis + // of the exact second derivatives below. + gInterpreter->ProcessLine(("clad::differentiate(" + namespaceName + "::roo_wrapper, \"input0[0]\");").c_str()); + + // The session tangent for the pushforward: the weights are constants, so + // their tangent is zero (a default-constructed Session holds the actual + // weights). Intermediate tensors inside this session are used as tangent + // scratch space during the forward pass. + _runtime->_zero_session.emplace(sessionName); + auto ptrZeroSession = toPtrString(_runtime->_zero_session.ptr, sessionName); + gInterpreter->ProcessLine((ptrZeroSession + "->SetWeightsToZero();").c_str()); + + // Directional-derivative wrapper: reverse differentiation of it w.r.t. + // both the inputs and the tangent direction gives the exact + // Hessian-vector product H . d_input (in the input adjoints) and the + // gradient (in the tangent adjoints). Both must be requested as active: + // for a non-varied argument clad would look up the custom helper pullbacks + // with reduced signatures, not find them, and silently fall back to + // differentiating the BLAS calls, which yields zero. + { + std::ostringstream ss; + ss << "namespace " << namespaceName << " {\n\n" + << "float roo_dir(Session const &session, Session const &zeroSession, " << innerParams << ", " + << innerTangentParams << ") {\n" + << " return roo_inner_wrapper_pushforward(session, " << innerArgs << ", zeroSession, " << innerTangentArgs + << ").pushforward;\n" + << "}\n\n" + << "float roo_dir_outer(Session const &session, Session const &zeroSession, " << innerParams << ", " + << innerTangentParams << ") {\n" + << " return roo_dir(session, zeroSession, " << innerArgs << ", " << innerTangentArgs << ");\n" + << "}\n\n" + << "} // namespace " << namespaceName << "\n"; + gInterpreter->Declare(ss.str().c_str()); + } + gInterpreter->ProcessLine( + ("clad::gradient(" + namespaceName + "::roo_dir_outer, \"" + cladInputs + ", " + innerTangentArgs + "\");") + .c_str()); + // The codegen call site (CodegenImpl::codegenImpl(RooONNXFunc)) passes one // double-array argument per input tensor. Emit roo_outer_wrapper and the matching // custom-derivative pullback with the corresponding number of parameters. @@ -422,7 +470,9 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: // as forward-mode derivatives that are then differentiated in reverse // mode: the forward pass needs a "pushforward" for the opaque // roo_outer_wrapper call, and the reverse pass over it needs the - // matching "pushforward_pullback". + // matching "pushforward_pullback". Both are exact: they are built on + // the clad-generated roo_inner_wrapper_pushforward and roo_dir_pullback + // (see above). // Comma-separated parameter helper for the emitted code. std::string tangentDoubleParams; // "double const *d_input0, ..." @@ -434,52 +484,42 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: tangentDoubleParams += "double const *d_input" + istr; } - // Emits "roo_outer_wrapper_pullback(0, ..., 1., 0, ...);". - auto emitGradCall = [&](std::string const &prefix, std::string const &gradPrefix) { - ss << " roo_outer_wrapper_pullback("; - for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << prefix << i << ", "; - } - ss << "1., "; - for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << gradPrefix << i << (i != nInputTensors - 1 ? ", " : ""); - } - ss << ");\n"; - }; - - // Emits zero-initialized double buffers named 0, ... with the - // per-tensor sizes. - auto emitBuffers = [&](std::string const &prefix) { + // Emits float conversion buffers for the inputs and the tangents. + auto emitFloatBuffers = [&]() { for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " double " << prefix << i << "[inputTensorDims[" << i << "].total_size()] = {};\n"; + ss << " float inputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " float dInputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n" + << " for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n" + << " inputFlt" << i << "[i] = input" << i << "[i];\n" + << " dInputFlt" << i << "[i] = d_input" << i << "[i];\n" + << " }\n"; } }; - // The pushforward evaluates the function value together with the - // directional derivative grad . d_input, both exactly via the custom - // pullback above. + // The pushforward evaluates the function value together with the exact + // directional derivative grad . d_input, in one forward pass. ss << "clad::ValueAndPushforward roo_outer_wrapper_pushforward(" << outerDoubleParams << ", " << tangentDoubleParams << ") {\n" - << " using namespace ::" << namespaceName << ";\n"; - emitBuffers("grad"); - emitGradCall("input", "grad"); - ss << " double dot = 0.;\n"; + << " using namespace ::" << namespaceName << ";\n" + << " auto &session = *" << ptrSession << ";\n" + << " auto &zeroSession = *" << ptrZeroSession << ";\n"; + emitFloatBuffers(); + ss << " auto vp = roo_inner_wrapper_pushforward(session, "; for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" - << " dot += grad" << i << "[i] * d_input" << i << "[i];\n" - << " }\n"; + ss << "inputFlt" << i << ", "; } - ss << " return {::" << namespaceName << "::roo_outer_wrapper(" << innerArgs << "), dot};\n" + ss << "zeroSession"; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << ", dInputFlt" << i; + } + ss << ");\n" + << " return {vp.value, vp.pushforward};\n" << "}\n\n"; // The pullback of the pushforward needs second derivatives only in the - // form of the Hessian-vector product H . d_input. SOFIE emits no - // forward-mode (pushforward) support for its operators, so instead of - // differentiating the model code again, the product is evaluated as a - // central finite difference of the *exact* generated gradient along the - // tangent direction. The step size (~cbrt of the float machine - // epsilon) balances the float-precision noise of the SOFIE gradient - // against the truncation error, giving ~1e-4 relative accuracy. + // form of the Hessian-vector product H . d_input, which roo_dir_pullback + // evaluates exactly (together with the gradient, as the adjoint of the + // tangent direction). ss << "void roo_outer_wrapper_pushforward_pullback(" << outerDoubleParams << ", " << tangentDoubleParams << ", clad::ValueAndPushforward d_y"; for (std::size_t i = 0; i < nInputTensors; ++i) { @@ -489,52 +529,41 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std:: ss << ", double *d_out_d_input" << i; } ss << ") {\n" - << " using namespace ::" << namespaceName << ";\n"; - emitBuffers("grad"); - emitGradCall("input", "grad"); - ss << " double scale = 1.;\n" - << " double norm = 0.;\n"; + << " using namespace ::" << namespaceName << ";\n" + << " auto &session = *" << ptrSession << ";\n" + << " auto &zeroSession = *" << ptrZeroSession << ";\n" + << " // Fresh adjoint sessions for every call: the clad-generated\n" + << " // second-order pullback does not restore the adjoint state inside\n" + << " // the sessions to zero (clad bug, see the \"clad referenced\n" + << " // '_tracker...' before its declaration\" warnings it prints during\n" + << " // generation), so reusing them would leak state between calls.\n" + << " Session dSession1;\n" + << " Session dSession2;\n"; + emitFloatBuffers(); for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" - << " scale = ::std::max(scale, ::std::abs(input" << i << "[i]));\n" - << " norm = ::std::max(norm, ::std::abs(d_input" << i << "[i]));\n" - << " }\n"; + ss << " float hvpFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n" + << " float gradFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n"; + } + ss << " roo_dir_pullback(session, zeroSession, "; + for (std::size_t i = 0; i < nInputTensors; ++i) { + ss << "inputFlt" << i << ", "; } - ss << " if (norm != 0.) {\n" - << " const double h = 6.7e-3 * scale;\n"; for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " double inputP" << i << "[inputTensorDims[" << i << "].total_size()];\n" - << " double inputM" << i << "[inputTensorDims[" << i << "].total_size()];\n" - << " for (::std::size_t i = 0; i < ::std::size(inputP" << i << "); ++i) {\n" - << " const double step = h * d_input" << i << "[i] / norm;\n" - << " inputP" << i << "[i] = input" << i << "[i] + step;\n" - << " inputM" << i << "[i] = input" << i << "[i] - step;\n" - << " }\n"; + ss << "dInputFlt" << i << ", "; } - ss << " "; - emitBuffers("gradP"); - ss << " "; - emitGradCall("inputP", "gradP"); - ss << " "; - emitBuffers("gradM"); - ss << " "; - emitGradCall("inputM", "gradM"); + ss << "1.F, &dSession1, &dSession2"; for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" - << " d_out_input" << i << "[i] += d_y.value * grad" << i << "[i]\n" - << " + d_y.pushforward * norm * (gradP" << i << "[i] - gradM" << i << "[i]) / (2. * h);\n" - << " }\n"; + ss << ", hvpFlt" << i; } - ss << " } else {\n"; for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" - << " d_out_input" << i << "[i] += d_y.value * grad" << i << "[i];\n" - << " }\n"; + ss << ", gradFlt" << i; } - ss << " }\n"; + ss << ");\n"; for (std::size_t i = 0; i < nInputTensors; ++i) { - ss << " for (::std::size_t i = 0; i < ::std::size(grad" << i << "); ++i) {\n" - << " d_out_d_input" << i << "[i] += d_y.pushforward * grad" << i << "[i];\n" + ss << " for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n" + << " d_out_input" << i << "[i] += d_y.value * gradFlt" << i << "[i] + d_y.pushforward * hvpFlt" << i + << "[i];\n" + << " d_out_d_input" << i << "[i] += d_y.pushforward * gradFlt" << i << "[i];\n" << " }\n"; } ss << "}\n\n" diff --git a/roofit/roofit/test/testRooONNXFunc.cxx b/roofit/roofit/test/testRooONNXFunc.cxx index a3c54838a8271..f2f2d5975261e 100644 --- a/roofit/roofit/test/testRooONNXFunc.cxx +++ b/roofit/roofit/test/testRooONNXFunc.cxx @@ -234,14 +234,13 @@ TEST(RooONNXFunc, CodegenHessian) roo_final.generateHessian(); - // The second derivatives are Hessian-vector products evaluated via finite - // differences of the exact gradient, which is itself limited by the float - // precision of the SOFIE computation. Hence the looser tolerance compared - // to the gradient checks. + // The second derivatives are exact Hessian-vector products (reverse-mode + // derivatives of the forward-mode derivative of the generated code), so + // they merit the same tolerance as the gradient checks. std::vector hess(n * n); roo_final.hessian(hess.data()); for (std::size_t i = 0; i < n * n; ++i) { - EXPECT_NEAR(hess[i], refHess[i], 1e-3); + EXPECT_NEAR(hess[i], refHess[i], 1e-5); } } @@ -273,7 +272,7 @@ TEST(RooONNXFunc, CodegenHessian_2Tensors) std::vector hess(nTotal * nTotal); roo_final.hessian(hess.data()); for (std::size_t i = 0; i < nTotal * nTotal; ++i) { - EXPECT_NEAR(hess[i], refHess[i], 1e-3); + EXPECT_NEAR(hess[i], refHess[i], 1e-5); } } @@ -310,7 +309,7 @@ TEST(RooONNXFunc, CodegenHessianCompound) for (std::size_t i = 0; i < n; ++i) { for (std::size_t j = 0; j < n; ++j) { const double ref = 2. * (refPred * refHess[i * n + j] + refGrad[i] * refGrad[j]); - EXPECT_NEAR(hess[i * n + j], ref, 1e-3); + EXPECT_NEAR(hess[i * n + j], ref, 1e-5); } } } From 2b1b21c9e0ca239afb1137af58936faa5bf92f5f Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 9 Oct 2026 09:03:21 +0000 Subject: [PATCH 7/9] [RF] Use Minuit2 directly in RooMinimizer, not via ROOT::Math::Minimizer MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RooMinimizer used to reach Minuit2 through the ROOT::Math::Minimizer plugin interface, while at the same time using Minuit2 directly as a library for things the generic interface can't express (vanishing second derivatives, the multiprocess gradient FCN). This commit removes that split: when the minimizer type is "Minuit2" (the default), RooMinimizer now uses the Minuit2 library interface directly (ModularFunctionMinimizer, MnHesse, MnMinos, MnContours, MnUserParameterState, FunctionMinimum). Other minimizer types still go through ROOT::Math::Minimizer. To make that possible, RooAbsMinimizerFcn now implements ROOT::Minuit2::FCNBase. RooMinimizerFcn provides gradient, Hessian and the vanishing-second-derivative mask through that interface, and the FCNBase adapter that MinuitFcnGrad carried for the multiprocess gradient is folded into MinuitFcnGrad itself. LikelihoodGradientJob gets the FCN from the RooMinimizer instead of from the Minuit2Minimizer wrapper. The internal FitResult is now the single source of truth for fit results on both paths; it additionally stores the global correlation coefficients, and exposes the covariance matrix status, which xRooFit uses instead of asking the (no longer existing) ROOT::Math::Minimizer. FitterInterface::GetMinimizer returns nullptr when Minuit2 is used directly. Behaviour is kept as in Minuit2Minimizer, including its status codes, print levels, strategy customization from the Minuit2 extra options, and the re-minimization when Minos finds a new minimum. Differences: Hesse now applies the current error definition, the "Minuit2Minimizer:" std::cout lines are now prefixed "RooMinimizer: Minuit2", and the Fumili algorithm, which never worked with RooFit FCNs, is rejected with an explicit error. 🤖 Done with the help of AI --- roofit/roofitcore/CMakeLists.txt | 4 +- roofit/roofitcore/inc/RooMinimizer.h | 25 +- roofit/roofitcore/src/RooAbsMinimizerFcn.h | 26 +- roofit/roofitcore/src/RooMinimizer.cxx | 835 +++++++++++++++++- roofit/roofitcore/src/RooMinimizerFcn.cxx | 50 +- roofit/roofitcore/src/RooMinimizerFcn.h | 13 +- .../TestStatistics/LikelihoodGradientJob.cxx | 14 +- .../src/TestStatistics/MinuitFcnGrad.cxx | 77 +- .../src/TestStatistics/MinuitFcnGrad.h | 18 +- roofit/roofitcore/test/testRooMinimizer.cxx | 2 +- roofit/xroofit/src/xRooFit.cxx | 12 +- 11 files changed, 934 insertions(+), 142 deletions(-) diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 4e2fe10fcd4b5..02527b41db31b 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -486,8 +486,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore ${EXTRA_DICT_OPTS} ) -# Minuit2 is needed in the implementation of RooMinimizerFcn, to advertise -# vanishing second derivatives to Minuit2Minimizer. +# RooMinimizer drives Minuit2 directly via its library interface, without +# going through the ROOT::Math::Minimizer plugin layer. target_link_libraries(RooFitCore PRIVATE Minuit2) # The following definitions are PUBLIC so they can also be used in ROOT-internal tests diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 4fae9157f263d..2238974c3bbfe 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -58,6 +58,7 @@ class RooMinimizer : public TObject { double Edm() const { return fEdm; } bool IsValid() const { return fValid; } int Status() const { return fStatus; } + int CovMatrixStatus() const { return fCovStatus; } void GetCovarianceMatrix(TMatrixDSym &cov) const; bool isParameterFixed(unsigned int ipar) const; @@ -120,6 +121,9 @@ class RooMinimizer : public TObject { } ROOT::Fit::FitConfig &Config() const { return *_config; } + /// The underlying ROOT::Math::Minimizer. Returns `nullptr` when the + /// minimizer type is Minuit2, because RooMinimizer then uses Minuit2 + /// directly without going through the ROOT::Math::Minimizer interface. ROOT::Math::Minimizer *GetMinimizer() const { return _minimizer; } const FitResult &Result() const { return *_result; } @@ -183,7 +187,8 @@ class RooMinimizer : public TObject { int evalCounter() const; void zeroEvalCount(); - /// Return underlying ROOT fitter object + /// Return an object giving access to the fit configuration and the + /// internal fit result, for backwards compatibility. inline auto fitter() { return std::make_unique(&_config, _minimizer.get(), _result.get()); } int getNPar() const; @@ -215,6 +220,7 @@ class RooMinimizer : public TObject { int exec(std::string const &algoName, std::string const &statusName); bool fitFCN(); + bool runMinimizer(); bool calculateHessErrors(); bool calculateMinosErrors(); @@ -229,10 +235,25 @@ class RooMinimizer : public TObject { void fillCorrMatrix(RooFitResult &fitRes); void updateErrors(); + // Direct interface to Minuit2 that bypasses ROOT::Math::Minimizer. It is + // used whenever the minimizer type is "Minuit2". + struct Minuit2State; + bool useMinuit2Directly() const { return _cfg.minimizerType == "Minuit2"; } + void initMinuit2(); + bool minuit2Minimize(); + bool minuit2Hesse(); + bool minuit2Minos(unsigned int index, double &errLow, double &errUp); + int runMinuit2Minos(unsigned int index, double &errLow, double &errUp); + bool + minuit2Contour(unsigned int ipar, unsigned int jpar, unsigned int npoints, double errorDef, double *x, double *y); + void setMinuit2ErrorDef(double up); + void fillResultFromMinuit2(bool isValid); + RooAbsReal &_function; ROOT::Fit::FitConfig _config; ///< fitter configuration (options and parameter settings) std::unique_ptr _result; /// _minimizer; /// _minimizer; /// _minuit2; /// +#include +#include #include #include #include #include // unique_ptr -class RooAbsMinimizerFcn { +/// Base class for the functions that RooMinimizer hands to Minuit. It +/// implements the ROOT::Minuit2::FCNBase interface, so that RooMinimizer can +/// drive Minuit2 directly. The function evaluation itself is implemented by +/// the derived classes. +class RooAbsMinimizerFcn : public ROOT::Minuit2::FCNBase { public: RooAbsMinimizerFcn(RooArgList paramList, RooMinimizer *context); - virtual ~RooAbsMinimizerFcn() = default; + ~RooAbsMinimizerFcn() override = default; + /// Set this function on a ROOT::Math::Minimizer. Only used for minimizers + /// other than Minuit2, which RooMinimizer uses directly via the FCNBase + /// interface that this class implements. virtual void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) = 0; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double Up() const override { return _errorDef; } + void SetErrorDef(double up) override { _errorDef = up; } + /// @} + /// Informs Minuit through its parameter_settings vector of RooFit parameter properties. bool synchronizeParameterSettings(std::vector ¶meters); /// Like synchronizeParameterSettings, Synchronize informs Minuit through @@ -122,6 +133,9 @@ class RooAbsMinimizerFcn { std::vector _floatableParamIndices; std::ofstream *_logfile = nullptr; + +private: + double _errorDef = 1.; ///< Minuit error definition, see FCNBase::Up(). }; #endif diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 133b70f7ffc77..9f32ff6ce674c 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -18,10 +18,11 @@ \class RooMinimizer \ingroup Roofitcore -Wrapper class around ROOT::Math::Minimizer that -provides a seamless interface between the minimizer functionality +RooMinimizer provides a seamless interface between the minimizer functionality and the native RooFit interface. -By default the Minimizer is Minuit 2. +By default the Minimizer is Minuit 2, which RooMinimizer drives directly via the +Minuit2 library interface. Other minimizer types are reached via the +ROOT::Math::Minimizer plugin interface. RooMinimizer can minimize any RooAbsReal function with respect to its parameters. Usual choices for minimization are the object returned by RooAbsPdf::createNLL() or RooAbsReal::createChi2(). @@ -64,11 +65,26 @@ Various methods are available to control verbosity or profiling. #include "RooFitImplHelpers.h" #include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include +#include +#include #include #include #include // logic_error @@ -138,8 +154,132 @@ std::unique_ptr setOperModesDirty(RooAbsReal &function) return {}; } +// Set the global Minuit2 print level for the duration of a Minuit2 call. +class Minuit2PrintLevelRAII { +public: + Minuit2PrintLevelRAII(int level) : _prevLevel{ROOT::Minuit2::MnPrint::SetGlobalLevel(level)} {} + ~Minuit2PrintLevelRAII() { ROOT::Minuit2::MnPrint::SetGlobalLevel(_prevLevel); } + +private: + int _prevLevel; +}; + +// The Minuit2-specific extra options, either set on the minimizer options +// explicitly or taken from the global defaults for Minuit2. +ROOT::Math::IOptions const *minuit2ExtraOptions(ROOT::Math::MinimizerOptions const &options) +{ + if (auto *opts = options.ExtraOptions()) { + return opts; + } + return ROOT::Math::MinimizerOptions::FindDefault("Minuit2"); +} + +// Build the Minuit2 strategy from the strategy level and the Minuit2-specific +// extra options. +ROOT::Minuit2::MnStrategy makeMinuit2Strategy(ROOT::Math::MinimizerOptions const &options) +{ + ROOT::Minuit2::MnStrategy st{static_cast(options.Strategy())}; + ROOT::Math::IOptions const *minuit2Opt = minuit2ExtraOptions(options); + if (!minuit2Opt) { + return st; + } + auto customize = [&minuit2Opt](const char *name, auto val) { + minuit2Opt->GetValue(name, val); + return val; + }; + st.SetGradientNCycles(customize("GradientNCycles", int(st.GradientNCycles()))); + st.SetHessianNCycles(customize("HessianNCycles", int(st.HessianNCycles()))); + st.SetHessianGradientNCycles(customize("HessianGradientNCycles", int(st.HessianGradientNCycles()))); + + st.SetGradientTolerance(customize("GradientTolerance", st.GradientTolerance())); + st.SetGradientStepTolerance(customize("GradientStepTolerance", st.GradientStepTolerance())); + st.SetHessianStepTolerance(customize("HessianStepTolerance", st.HessianStepTolerance())); + st.SetHessianG2Tolerance(customize("HessianG2Tolerance", st.HessianG2Tolerance())); + + st.SetHessianCentralFDMixedDerivatives( + customize("HessianCentralFDMixedDerivatives", int(st.HessianCentralFDMixedDerivatives()))); + st.SetHessianForcePosDef(customize("HessianForcePosDef", int(st.HessianForcePosDef()))); + + return st; +} + +// Create the Minuit2 minimization algorithm for the given algorithm name. +// Returns nullptr for algorithms that RooFit can't use. +std::unique_ptr makeMinuit2Minimizer(std::string algo) +{ + std::transform(algo.begin(), algo.end(), algo.begin(), [](unsigned char c) { return std::tolower(c); }); + using namespace ROOT::Minuit2; + if (algo == "simplex") + return std::make_unique(); + if (algo == "minimize") + return std::make_unique(); + if (algo == "scan") + return std::make_unique(); + if (algo == "bfgs") + return std::make_unique(VariableMetricMinimizer::BFGSType()); + if (algo == "fumili" || algo == "fumili2") { + // Fumili needs a FumiliFCNBase, which RooFit does not provide. + return nullptr; + } + // MIGRAD is the default algorithm, also for unknown algorithm names. + return std::make_unique(); +} + +bool isFixedOrConst(ROOT::Minuit2::MinuitParameter const &par) +{ + return par.IsFixed() || par.IsConst(); +} + +// Covariance matrix element in the external parameter space, zero for fixed +// parameters and when no covariance matrix is available. +double minuit2Covariance(ROOT::Minuit2::MnUserParameterState const &state, unsigned int i, unsigned int j) +{ + if (!state.HasCovariance() || isFixedOrConst(state.Parameter(i)) || isFixedOrConst(state.Parameter(j))) { + return 0.; + } + return state.Covariance()(state.IntOfExt(i), state.IntOfExt(j)); +} + +// Covariance matrix status code with the convention of ROOT::Math::Minimizer: +// -1: not available (inversion failed or Hesse failed) +// 0: available but not positive defined +// 1: covariance only approximate +// 2: full matrix but forced positive definite +// 3: full accurate matrix +int minuit2CovMatrixStatus(ROOT::Minuit2::FunctionMinimum const *minimum, + ROOT::Minuit2::MnUserParameterState const &state) +{ + if (!minimum) { + return state.CovarianceStatus(); + } + if (minimum->HasAccurateCovar()) + return 3; + if (minimum->HasMadePosDefCovar()) + return 2; + if (minimum->HasValidCovariance()) + return 1; + if (minimum->HasCovariance()) + return 0; + return -1; +} + } // namespace +/// State of the direct Minuit2 interface. +struct RooMinimizer::Minuit2State { + /// Minuit2 parameter state: the starting point for the next minimization, + /// and the result of the last operation. + ROOT::Minuit2::MnUserParameterState state; + /// Result of the last minimization, needed for Hesse, Minos and contours. + std::unique_ptr minimum; + /// Status code with the same convention as ROOT::Minuit2::Minuit2Minimizer. + int status = 0; + /// Status of the last Minos run. + int minosStatus = -1; + /// Whether the parameter state was initialized from the parameter settings. + bool initialized = false; +}; + //////////////////////////////////////////////////////////////////////////////// /// Construct MINUIT interface to given function. Function can be anything, /// but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2 @@ -211,6 +351,7 @@ void RooMinimizer::initMinimizerFirstPart() { RooSentinel::activate(); setMinimizerType(""); + _minuit2 = std::make_unique(); _config.SetMinimizer(_cfg.minimizerType.c_str()); setEps(1.0); // default tolerance @@ -428,7 +569,7 @@ int RooMinimizer::exec(std::string const &algoName, std::string const &statusNam int RooMinimizer::hesse() { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::hesse: Error, run Migrad before Hesse!" << std::endl; _status = -1; return _status; @@ -445,7 +586,7 @@ int RooMinimizer::hesse() int RooMinimizer::minos() { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!" << std::endl; _status = -1; return _status; @@ -462,7 +603,7 @@ int RooMinimizer::minos() int RooMinimizer::minos(const RooArgSet &minosParamList) { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!" << std::endl; _status = -1; } else if (!minosParamList.empty()) { @@ -570,7 +711,7 @@ void RooMinimizer::optimizeConst(int /*flag*/) {} RooFit::OwningPtr RooMinimizer::save(const char *userName, const char *userTitle) { - if (_minimizer == nullptr) { + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::save: Error, run minimization before!" << std::endl; return nullptr; } @@ -584,7 +725,7 @@ RooFit::OwningPtr RooMinimizer::save(const char *userName, const c fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()); fitRes->setStatus(_status); - fitRes->setCovQual(_minimizer->CovMatrixStatus()); + fitRes->setCovQual(_result->fCovStatus); fitRes->setMinNLL(_result->fVal - _fcn->getOffset()); fitRes->setEDM(_result->fEdm); @@ -628,7 +769,7 @@ void RooMinimizer::fillCorrMatrix(RooFitResult &fitRes) const std::size_t nParams = _fcn->getNDim(); TMatrixDSym corrs(nParams); TMatrixDSym covs(nParams); - std::vector globalCC = _minimizer->GlobalCC(); + std::vector globalCC = _result->fGlobalCC; globalCC.resize(nParams); // pad with zeros for (std::size_t ic = 0; ic < nParams; ic++) { for (std::size_t ii = 0; ii < nParams; ii++) { @@ -684,15 +825,23 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do TMarker *point = new TMarker(var1.getVal(), var2.getVal(), 8); frame->addObject(point); - // check first if a inimizer is available. If not means - // the minimization is not done , so do it - if (_minimizer == nullptr) { + // check first if a minimization was done + if (_result == nullptr) { coutW(Minimization) << "RooMinimizer::contour: Error, run Migrad before contours!" << std::endl; return frame; } // remember our original value of ERRDEF - double errdef = _minimizer->ErrorDef(); + const double errdef = _config.MinimizerOptions().ErrorDef(); + + // compute a contour at the given error level + auto computeContour = [&](double up, double *xcoor, double *ycoor) { + if (useMinuit2Directly()) { + return minuit2Contour(index1, index2, npoints, up, xcoor, ycoor); + } + _minimizer->SetErrorDef(up); + return _minimizer->Contour(index1, index2, npoints, xcoor, ycoor); + }; double n[6]; n[0] = n1; @@ -706,13 +855,10 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do for (int ic = 0; ic < 6; ic++) { if (n[ic] > 0) { - // set the value corresponding to an n1-sigma contour - _minimizer->SetErrorDef(n[ic] * n[ic] * errdef); - - // calculate and draw the contour + // calculate and draw the contour at the level corresponding to n-sigma std::vector xcoor(npoints + 1); std::vector ycoor(npoints + 1); - bool ret = _minimizer->Contour(index1, index2, npoints, xcoor.data(), ycoor.data()); + bool ret = computeContour(n[ic] * n[ic] * errdef, xcoor.data(), ycoor.data()); if (!ret) { coutE(Minimization) << "RooMinimizer::contour(" << GetName() @@ -734,7 +880,11 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do } // restore the original ERRDEF - _minimizer->SetErrorDef(errdef); + if (useMinuit2Directly()) { + setMinuit2ErrorDef(errdef); + } else { + _minimizer->SetErrorDef(errdef); + } // restore parameter values params.assign(paramSave); @@ -805,7 +955,7 @@ RooFit::OwningPtr RooMinimizer::lastMinuitFit() // Import the results of the last fit performed, interpreting // the fit parameters as the given varList of parameters. - if (_minimizer == nullptr) { + if (_result == nullptr) { oocoutE(nullptr, InputArguments) << "RooMinimizer::save: Error, run minimization before!" << std::endl; return nullptr; } @@ -851,7 +1001,7 @@ RooFit::OwningPtr RooMinimizer::lastMinuitFit() res->setFinalParList(floatPars); res->setMinNLL(_result->fVal); res->setEDM(_result->fEdm); - res->setCovQual(_minimizer->CovMatrixStatus()); + res->setCovQual(_result->fCovStatus); res->setStatus(_result->fStatus); fillCorrMatrix(*res); @@ -954,10 +1104,7 @@ bool RooMinimizer::fitFCN() if (nPdfs == 0) { coutI(Minimization) << "[fitFCN] No discrete parameters, performing continuous minimization only" << std::endl; FreezeDisconnectedParametersRAII freeze(this, *_fcn); - bool isValid = _minimizer->Minimize(); - if (!_result) - _result = std::make_unique(); - fillResult(isValid); + bool isValid = runMinimizer(); if (isValid) updateFitConfig(); return isValid; @@ -994,12 +1141,12 @@ bool RooMinimizer::fitFCN() pdfIndices[i]->setConstant(true); } FreezeDisconnectedParametersRAII freeze(this, *_fcn); - _minimizer->Minimize(); + runMinimizer(); for (size_t i = 0; i < nPdfs; ++i) pdfIndices[i]->setConstant(wasConst[i]); - double val = _minimizer->MinValue(); + double val = _result->fVal; tried.insert(combo); nllMap[combo] = val; @@ -1016,7 +1163,7 @@ bool RooMinimizer::fitFCN() } FreezeDisconnectedParametersRAII freeze(this, *_fcn); - _minimizer->Minimize(); + runMinimizer(); coutI(Minimization) << "All NLL Values per Combination:" << std::endl; for (const auto &entry : nllMap) { @@ -1046,13 +1193,24 @@ bool RooMinimizer::fitFCN() coutI(Minimization) << ssBest.str() << std::endl; - if (!_result) - _result = std::make_unique(); - fillResult(true); + _result->fValid = true; updateFitConfig(); return true; } + +/// Run the minimization algorithm from the current minimizer state and store +/// the outcome in the internal fit result. Returns whether a valid minimum was +/// found. +bool RooMinimizer::runMinimizer() +{ + bool isValid = useMinuit2Directly() ? minuit2Minimize() : _minimizer->Minimize(); + if (!_result) + _result = std::make_unique(); + fillResult(isValid); + return isValid; +} + bool RooMinimizer::calculateHessErrors() { // compute the Hesse errors according to configuration @@ -1060,14 +1218,17 @@ bool RooMinimizer::calculateHessErrors() auto operModeRAII = setOperModesDirty(_function); - // update minimizer (recreate if not done or if name has changed - if (!updateMinimizerOptions()) { - coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; - return false; + bool ret = false; + if (useMinuit2Directly()) { + ret = minuit2Hesse(); + } else { + // update minimizer (recreate if not done or if name has changed + if (!updateMinimizerOptions()) { + coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; + return false; + } + ret = _minimizer->Hesse(); } - - // run Hesse - bool ret = _minimizer->Hesse(); if (!ret) coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error when calculating Hessian" << std::endl; @@ -1097,8 +1258,8 @@ bool RooMinimizer::calculateMinosErrors() auto operModeRAII = setOperModesDirty(_function); // update minimizer (but cannot re-create in this case). Must use an existing one - if (!updateMinimizerOptions(false)) { - coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; + if (!useMinuit2Directly() && !updateMinimizerOptions(false)) { + coutE(Minimization) << "RooMinimizer::calculateMinosErrors() Error re-initializing the minimizer" << std::endl; return false; } @@ -1119,9 +1280,17 @@ bool RooMinimizer::calculateMinosErrors() for (int i = 0; i < iparMax; ++i) { double elow, eup; unsigned int index = (!ipars.empty()) ? ipars[i] : i; - bool ret = _minimizer->GetMinosError(index, elow, eup); + bool ret = false; + int minosStatus = 0; + if (useMinuit2Directly()) { + ret = minuit2Minos(index, elow, eup); + minosStatus = _minuit2->minosStatus; + } else { + ret = _minimizer->GetMinosError(index, elow, eup); + minosStatus = _minimizer->MinosStatus(); + } // flags case when a new minimum has been found - if ((_minimizer->MinosStatus() & 8) != 0) { + if ((minosStatus & 8) != 0) { iparNewMin = i; } if (ret) @@ -1147,6 +1316,10 @@ bool RooMinimizer::calculateMinosErrors() void RooMinimizer::initMinimizer() { + if (useMinuit2Directly()) { + initMinuit2(); + return; + } _minimizer = std::unique_ptr(_config.CreateMinimizer()); _fcn->initMinimizer(*_minimizer, this); _minimizer->SetVariables(_config.ParamsSettings().begin(), _config.ParamsSettings().end()); @@ -1233,6 +1406,10 @@ RooMinimizer::FitResult::FitResult(const ROOT::Fit::FitConfig &fconfig) void RooMinimizer::fillResult(bool isValid) { + if (useMinuit2Directly()) { + fillResultFromMinuit2(isValid); + return; + } ROOT::Math::Minimizer &min = *_minimizer; ROOT::Fit::FitConfig const &fconfig = _config; @@ -1270,6 +1447,7 @@ void RooMinimizer::fillResult(bool isValid) // if minimizer provides error provides also error matrix // clear in case of re-filling an existing result _result->fCovMatrix.clear(); + _result->fGlobalCC.clear(); if (min.Errors() != nullptr) { updateErrors(); @@ -1278,6 +1456,10 @@ void RooMinimizer::fillResult(bool isValid) bool RooMinimizer::update(bool isValid) { + if (useMinuit2Directly()) { + fillResultFromMinuit2(isValid); + return true; + } ROOT::Math::Minimizer &min = *_minimizer; ROOT::Fit::FitConfig const &fconfig = _config; @@ -1312,6 +1494,7 @@ void RooMinimizer::updateErrors() _result->fErrors.resize(npar); std::copy(min.Errors(), min.Errors() + npar, _result->fErrors.begin()); + _result->fGlobalCC = min.GlobalCC(); if (_result->fCovStatus != 0) { @@ -1327,6 +1510,574 @@ void RooMinimizer::updateErrors() // minos errors are set separately when calling Fitter::CalculateMinosErrors() } +//////////////////////////////////////////////////////////////////////////////// +/// Initialize the Minuit2 parameter state from the parameter settings, and +/// forget about any previous function minimum. + +void RooMinimizer::initMinuit2() +{ + using ROOT::Minuit2::MnUserParameterState; + + _minuit2->state = MnUserParameterState{}; + _minuit2->minimum.reset(); + _minuit2->status = 0; + _minuit2->minosStatus = -1; + + MnUserParameterState &state = _minuit2->state; + + for (ROOT::Fit::ParameterSettings const &par : _config.ParamsSettings()) { + if (par.IsFixed()) { + // Fixed parameters still need a step size, otherwise Minuit2 would + // treat them as constants that can't be released anymore. + const double step = par.Value() != 0 ? 0.1 * std::abs(par.Value()) : 0.1; + state.Add(par.Name(), par.Value(), step); + state.Fix(par.Name()); + continue; + } + if (par.StepSize() <= 0) { + // A parameter without a valid step size is treated as a constant. + state.Add(par.Name(), par.Value()); + } else { + state.Add(par.Name(), par.Value(), par.StepSize()); + } + if (par.IsDoubleBound()) { + state.SetLimits(par.Name(), par.LowerLimit(), par.UpperLimit()); + } else if (par.HasLowerLimit()) { + state.SetLowerLimit(par.Name(), par.LowerLimit()); + } else if (par.HasUpperLimit()) { + state.SetUpperLimit(par.Name(), par.UpperLimit()); + } + } + + if (_cfg.setInitialCovariance) { + // Diagonal covariance matrix from the parameter errors, in the packed + // lower-triangular storage of MnUserCovariance. + const unsigned int n = _fcn->getNDim(); + ROOT::Minuit2::MnUserCovariance cov{n}; + for (unsigned int i = 0; i < n; ++i) { + const double err = _fcn->floatableParam(i).getError(); + cov(i, i) = err * err; + } + state.AddCovariance(cov); + } + + _minuit2->initialized = true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Set the error definition on the function and on the function minimum, if +/// there is one. + +void RooMinimizer::setMinuit2ErrorDef(double up) +{ + _fcn->SetErrorDef(up); + if (_minuit2->minimum && _minuit2->minimum->Up() != up) { + _minuit2->minimum->SetErrorDef(up); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run the configured Minuit2 minimization algorithm, starting from the current +/// parameter state. Returns whether a valid minimum was found. + +bool RooMinimizer::minuit2Minimize() +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + + if (!_minuit2->initialized) { + initMinuit2(); + } + + auto minimizer = makeMinuit2Minimizer(options.MinimizerAlgorithm()); + if (!minimizer) { + coutE(Minimization) << "RooMinimizer::minimize: the Minuit2 algorithm \"" << options.MinimizerAlgorithm() + << "\" is not supported by RooFit" << std::endl; + return false; + } + + MnUserParameterState &state = _minuit2->state; + + // delete result of previous minimization + _minuit2->minimum.reset(); + + const unsigned int maxfcn = options.MaxFunctionCalls(); + const double tol = options.Tolerance(); + const int printLevel = options.PrintLevel(); + _fcn->SetErrorDef(options.ErrorDef()); + + if (printLevel >= 1) { + // print the real number of maxfcn used (defined in ModularFunctionMinimizer) + unsigned int maxfcnUsed = maxfcn; + if (maxfcnUsed == 0) { + const unsigned int nvar = state.VariableParameters(); + maxfcnUsed = 200 + 100 * nvar + 5 * nvar * nvar; + } + std::cout << "RooMinimizer: Minuit2 minimize with max-calls " << maxfcnUsed << " convergence for edm < " << tol + << " strategy " << options.Strategy() << std::endl; + } + + minimizer->Builder().SetPrintLevel(printLevel); + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + if (ROOT::Math::IOptions const *minuit2Opt = minuit2ExtraOptions(options)) { + int storageLevel = 1; + if (minuit2Opt->GetValue("StorageLevel", storageLevel)) { + minimizer->Builder().SetStorageLevel(storageLevel); + } + if (printLevel > 0) { + std::cout << "RooMinimizer: Minuit2 - Changing default options" << std::endl; + minuit2Opt->Print(); + } + } + + const MnStrategy strategy = makeMinuit2Strategy(options); + + _minuit2->minimum = std::make_unique(minimizer->Minimize(*_fcn, state, strategy, maxfcn, tol)); + FunctionMinimum const &minimum = *_minuit2->minimum; + + // copy minimum state (parameter values and errors) + state = minimum.UserState(); + + // Determine the status code, with the same convention as Minuit2Minimizer. + int &status = _minuit2->status; + status = 0; + std::string txt; + if (!minimum.HasPosDefCovar()) { + // this happens normally when Hesse failed + // it can happen in case MnSeed failed (see ROOT-9522) + txt = "Covar is not pos def"; + status = 5; + } + if (minimum.HasMadePosDefCovar()) { + txt = "Covar was made pos def"; + status = 1; + } + if (minimum.HesseFailed()) { + txt = "Hesse is not valid"; + status = 2; + } + if (minimum.IsAboveMaxEdm()) { + txt = "Edm is above max"; + status = 3; + } + if (minimum.HasReachedCallLimit()) { + txt = "Reached call limit"; + status = 4; + } + + MnPrint print("RooMinimizer::minimize", printLevel); + const bool validMinimum = minimum.IsValid(); + if (validMinimum) { + // print a warning message in case something is not ok + if (status != 0 && printLevel > 0) + print.Warn(txt); + } else { + // minimum is not valid when state is not valid and edm is over max or has passed call limits + if (status == 0) { + // this should not happen + txt = "unknown failure"; + status = 6; + } + print.Warn("Minimization did NOT converge,", txt); + } + + if (printLevel >= 1) { + std::cout << "RooMinimizer: Minuit2 " << (validMinimum ? "valid" : "invalid") << " minimum - status = " << status + << std::endl; + const int prec = std::cout.precision(18); + std::cout << "FVAL = " << state.Fval() << std::endl; + std::cout << "Edm = " << state.Edm() << std::endl; + std::cout.precision(prec); + std::cout << "Nfcn = " << state.NFcn() << std::endl; + if (validMinimum) { + for (MinuitParameter const &par : state.MinuitParameters()) { + std::cout << par.Name() << "\t = " << par.Value() << "\t "; + if (par.IsFixed()) + std::cout << "(fixed)" << std::endl; + else if (par.IsConst()) + std::cout << "(const)" << std::endl; + else if (par.HasLimits()) + std::cout << "+/- " << par.Error() << "\t(limited)" << std::endl; + else + std::cout << "+/- " << par.Error() << std::endl; + } + } + } + + return validMinimum; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run Minuit2's HESSE. If a function minimum from a previous minimization +/// exists, it is updated with the result. Returns whether a valid covariance +/// matrix was obtained. + +bool RooMinimizer::minuit2Hesse() +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + + if (!_minuit2->initialized) { + initMinuit2(); + } + + const unsigned int maxfcn = options.MaxFunctionCalls(); + const int printLevel = options.PrintLevel(); + + MnPrint print("RooMinimizer::hesse", printLevel); + print.Info("Using max-calls", maxfcn); + + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + MnUserParameterState &state = _minuit2->state; + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + setMinuit2ErrorDef(options.ErrorDef()); + + MnHesse hesse(makeMinuit2Strategy(options)); + + if (_minuit2->minimum) { + // run hesse and function minimum will be updated with Hesse result + hesse(*_fcn, *_minuit2->minimum, maxfcn); + state = _minuit2->minimum->UserState(); + } else { + // run Hesse on point stored in current state (independent of function minimum validity) + state = hesse(*_fcn, state, maxfcn); + } + + if (printLevel >= 3) { + std::cout << "RooMinimizer::hesse - State returned from Hesse " << std::endl; + std::cout << state << std::endl; + } + + const int covStatus = state.CovarianceStatus(); + std::string covStatusType = "not valid"; + if (covStatus == 1) + covStatusType = "approximate"; + if (covStatus == 2) + covStatusType = "full but made positive defined"; + if (covStatus == 3) + covStatusType = "accurate"; + if (covStatus == 0) + covStatusType = "full but not positive defined"; + + if (!state.HasCovariance()) { + // if false means error is not valid and this is due to a failure in Hesse + // update minimizer error status + int hstatus = 4; + // information on error state can be retrieved only if the minimum is available + if (_minuit2->minimum) { + if (_minuit2->minimum->Error().HesseFailed()) + hstatus = 1; + if (_minuit2->minimum->Error().InvertFailed()) + hstatus = 2; + else if (!(_minuit2->minimum->Error().IsPosDef())) + hstatus = 3; + } + + print.Warn("Hesse failed - matrix is", covStatusType); + print.Warn(hstatus); + + _minuit2->status += 100 * hstatus; + return false; + } + + print.Info("Hesse is valid - matrix is", covStatusType); + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run Minuit2's MINOS for the parameter with the given index. If a new +/// minimum is found, the minimization is repeated from there and MINOS is run +/// again. Returns whether both errors are valid. + +bool RooMinimizer::minuit2Minos(unsigned int index, double &errLow, double &errUp) +{ + using namespace ROOT::Minuit2; + + errLow = 0; + errUp = 0; + + const int printLevel = _config.MinimizerOptions().PrintLevel(); + MnPrint print("RooMinimizer::minos", printLevel); + + if (!_minuit2->minimum) { + print.Error("Failed - no function minimum existing"); + return false; + } + + MnUserParameterState &state = _minuit2->state; + + // need to know if parameter is const or fixed + if (isFixedOrConst(state.Parameter(index))) { + return false; + } + + if (!_minuit2->minimum->IsValid()) { + print.Error("Failed - invalid function minimum"); + return false; + } + + setMinuit2ErrorDef(_config.MinimizerOptions().ErrorDef()); + + int mstatus = runMinuit2Minos(index, errLow, errUp); + + // run again the Minimization in case of a new minimum + // bit 8 is set + if ((mstatus & 8) != 0) { + print.Info([&](std::ostream &os) { + os << "Found a new minimum: run again the Minimization starting from the new point"; + os << "\nFVAL = " << state.Fval(); + for (MinuitParameter const &par : state.MinuitParameters()) { + os << '\n' << par.Name() << "\t = " << par.Value(); + } + }); + // release parameter that was fixed in the returned state from Minos + state.Release(index); + if (!minuit2Minimize()) + return false; + // run again Minos from new Minimum (also lower error needs to be re-computed) + print.Info("Run now again Minos from the new found Minimum"); + mstatus = runMinuit2Minos(index, errLow, errUp); + + // do not reset new minimum bit to flag for other parameters + mstatus |= 8; + } + + _minuit2->status += 10 * mstatus; + _minuit2->minosStatus = mstatus; + + return ((mstatus & 1) == 0) && ((mstatus & 2) == 0); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Run MINOS once for the parameter with the given index and return the MINOS +/// status bits, with the same convention as Minuit2Minimizer: +/// bit 1: lower error invalid, bit 2: upper error invalid, bit 4: invalid +/// because the maximum number of function calls was reached, bit 8: invalid +/// because a new minimum was found, bit 16: parameter is at a limit. + +int RooMinimizer::runMinuit2Minos(unsigned int index, double &errLow, double &errUp) +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + const int printLevel = options.PrintLevel(); + MnPrint print("RooMinimizer::minos", printLevel); + + Minuit2PrintLevelRAII printLevelRAII{printLevel}; + + MnUserParameterState &state = _minuit2->state; + if (options.Precision() > 0) { + state.SetPrecision(options.Precision()); + } + + // Like Minuit2Minimizer, use the default strategy for Minos. + MnMinos minos(*_fcn, *_minuit2->minimum); + + const unsigned int maxfcn = options.MaxFunctionCalls(); + // Tolerance for the migrad calls inside Minos. Cut off too small values, + // which are not needed. + const double tol = std::max(options.Tolerance(), 0.01); + + const char *parName = state.Name(index); + + if (printLevel >= 1) { + // get the real number of maxfcn used (defined in MnMinos) to be printed + unsigned int maxfcnUsed = maxfcn; + if (maxfcnUsed == 0) { + const unsigned int nvar = state.VariableParameters(); + maxfcnUsed = 2 * (nvar + 1) * (200 + 100 * nvar + 5 * nvar * nvar); + } + std::cout << "RooMinimizer::minos - Run MINOS for parameter #" << index << " : " << parName << " using max-calls " + << maxfcnUsed << ", tolerance " << tol << std::endl; + } + + const MinosError me = minos.Minos(index, maxfcn, tol); + + // Note that the only invalid condition can happen when the (npar-1) minimization fails + // The error is also invalid when the maximum number of calls is reached or a new function minimum is found + // in case of the parameter at the limit the error is not invalid. + // When the error is invalid the returned error is the Hessian error. + if (!me.LowerValid()) { + print.Warn("Invalid lower error for parameter", parName); + } + if (!me.UpperValid()) { + print.Warn("Invalid upper error for parameter", parName); + } + if (me.AtLowerLimit()) { + print.Warn("Lower error for parameter", parName, "is at the Lower limit!"); + } + if (me.AtUpperLimit()) { + print.Warn("Upper error for parameter", parName, "is at the Upper limit!"); + } + if (me.AtLowerMaxFcn()) { + print.Warn("Maximum number of function calls exceeded when running for lower error for parameter", parName); + } + if (me.AtUpperMaxFcn()) { + print.Warn("Maximum number of function calls exceeded when running for upper error for parameter", parName); + } + if (me.LowerNewMin()) { + print.Warn("New Minimum found while running Minos for lower error for parameter", parName); + } + if (me.UpperNewMin()) { + print.Warn("New Minimum found while running Minos for upper error for parameter", parName); + } + if (printLevel >= 1) { + if (me.LowerValid()) + std::cout << "Minos: Lower error for parameter " << parName << " : " << me.Lower() << std::endl; + if (me.UpperValid()) + std::cout << "Minos: Upper error for parameter " << parName << " : " << me.Upper() << std::endl; + } + + int mstatus = 0; + if (!me.LowerValid()) { + mstatus |= 1; + if (me.AtLowerMaxFcn()) + mstatus |= 4; + if (me.LowerNewMin()) + mstatus |= 8; + } + if (!me.UpperValid()) { + mstatus |= 2; + if (me.AtUpperMaxFcn()) + mstatus |= 4; + if (me.UpperNewMin()) + mstatus |= 8; + } + if (me.AtUpperLimit() || me.AtLowerLimit()) + mstatus |= 16; + + errLow = me.Lower(); + errUp = me.Upper(); + + // in case of new minimum found update also the minimum state + if (me.LowerNewMin() && me.UpperNewMin()) { + // take state with lower function value + state = (me.LowerState().Fval() < me.UpperState().Fval()) ? me.LowerState() : me.UpperState(); + } else if (me.LowerNewMin()) { + state = me.LowerState(); + } else if (me.UpperNewMin()) { + state = me.UpperState(); + } + + return mstatus; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Compute a contour for parameters ipar and jpar at the given error level +/// with Minuit2's MnContours, which requires a valid function minimum. + +bool RooMinimizer::minuit2Contour(unsigned int ipar, unsigned int jpar, unsigned int npoints, double errorDef, + double *x, double *y) +{ + using namespace ROOT::Minuit2; + + ROOT::Math::MinimizerOptions const &options = _config.MinimizerOptions(); + const int printLevel = options.PrintLevel(); + MnPrint print("RooMinimizer::contour", printLevel); + + if (!_minuit2->minimum) { + print.Error("No function minimum existing; must minimize function before"); + return false; + } + + if (!_minuit2->minimum->IsValid()) { + print.Error("Invalid function minimum"); + return false; + } + + setMinuit2ErrorDef(errorDef); + + print.Info("Computing contours at level -", errorDef); + + // switch off Minuit2 printing (for level of 0,1) + Minuit2PrintLevelRAII printLevelRAII{printLevel - 1}; + + if (options.Precision() > 0) { + _minuit2->state.SetPrecision(options.Precision()); + } + + // eventually one should specify tolerance in contours + MnContours contour(*_fcn, *_minuit2->minimum, makeMinuit2Strategy(options)); + + std::vector> result = contour(ipar, jpar, npoints); + if (result.size() != npoints) { + print.Error("Invalid result from MnContours"); + return false; + } + for (unsigned int i = 0; i < npoints; ++i) { + x[i] = result[i].first; + y[i] = result[i].second; + } + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fill the internal fit result from the Minuit2 parameter state. + +void RooMinimizer::fillResultFromMinuit2(bool isValid) +{ + using namespace ROOT::Minuit2; + + MnUserParameterState const &state = _minuit2->state; + std::vector const ¶ms = state.MinuitParameters(); + const unsigned int npar = params.size(); + + _result->fValid = isValid; + _result->fStatus = _minuit2->status; + _result->fCovStatus = minuit2CovMatrixStatus(_minuit2->minimum.get(), state); + _result->fVal = state.Fval(); + _result->fEdm = state.Edm(); + _result->fMinimType = _config.MinimizerName(); + + _result->fParams.resize(npar); + _result->fErrors.resize(npar); + for (unsigned int i = 0; i < npar; ++i) { + _result->fParams[i] = params[i].Value(); + _result->fErrors[i] = isFixedOrConst(params[i]) ? 0. : params[i].Error(); + } + + // check for fixed parameters + for (unsigned int ipar = 0; ipar < npar; ++ipar) { + if (_config.ParSettings(ipar).IsFixed()) + _result->fFixedParams[ipar] = true; + } + + // fill error matrix in packed lower-triangular storage + _result->fCovMatrix.clear(); + if (_result->fCovStatus != 0) { + _result->fCovMatrix.reserve(npar * (npar + 1) / 2); + for (unsigned int i = 0; i < npar; ++i) { + for (unsigned int j = 0; j <= i; ++j) { + _result->fCovMatrix.push_back(minuit2Covariance(state, i, j)); + } + } + } + + // global correlation coefficients, zero for fixed parameters + _result->fGlobalCC.clear(); + MnGlobalCorrelationCoeff const globalCC = state.GlobalCC(); + if (globalCC.IsValid()) { + _result->fGlobalCC.resize(npar); + for (unsigned int i = 0; i < npar; ++i) { + _result->fGlobalCC[i] = isFixedOrConst(params[i]) ? 0. : globalCC.GlobalCC()[state.IntOfExt(i)]; + } + } + // minos errors are set separately when calling calculateMinosErrors() +} + double RooMinimizer::FitResult::lowerError(unsigned int i) const { // return lower Minos error for parameter i diff --git a/roofit/roofitcore/src/RooMinimizerFcn.cxx b/roofit/roofitcore/src/RooMinimizerFcn.cxx index 6a3235950f2b3..efe06b57e8184 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooMinimizerFcn.cxx @@ -16,8 +16,10 @@ ////////////////////////////////////////////////////////////////////////////// /// \class RooMinimizerFcn -/// RooMinimizerFcn is an interface to the ROOT::Math::IBaseFunctionMultiDim, -/// a function that ROOT's minimisers use to carry out minimisations. +/// RooMinimizerFcn adapts a RooAbsReal to the function interface of the +/// minimizers used by RooMinimizer: the ROOT::Minuit2::FCNBase interface for +/// Minuit2, and ROOT::Math::IBaseFunctionMultiDim for the other minimizers +/// that are reached via ROOT::Math::Minimizer. /// #include "RooMinimizerFcn.h" @@ -35,7 +37,6 @@ #include "RooRealVar.h" #include "Math/Functor.h" -#include "Minuit2/Minuit2Minimizer.h" #include "TMatrixDSym.h" #include @@ -151,14 +152,18 @@ RooMinimizerFcn::RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context) { unsigned int nDim = getNDim(); - if (context->_cfg.useGradient && funct->hasGradient()) { + _useGradient = context->_cfg.useGradient && funct->hasGradient(); + _useHessian = context->_cfg.useHessian && funct->hasHessian(); + + auto evalFunc = [this](const double *x) { return (*this)(x); }; + if (_useGradient) { _gradientOutput.resize(_allParams.size()); - _multiGenFcn = std::make_unique(this, &RooMinimizerFcn::operator(), - &RooMinimizerFcn::evaluateGradient, nDim); + auto gradFunc = [this](const double *x, double *out) { evaluateGradient(x, out); }; + _multiGenFcn = std::make_unique(evalFunc, nDim, gradFunc); } else { - _multiGenFcn = std::make_unique(std::cref(*this), nDim); + _multiGenFcn = std::make_unique(evalFunc, nDim); } - if (context->_cfg.useHessian) { + if (_useHessian) { _hessianOutput.resize(_allParams.size() * _allParams.size()); } } @@ -311,21 +316,28 @@ bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) co return true; } -void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer *context) +std::vector RooMinimizerFcn::Gradient(std::vector const &x) const +{ + std::vector out(getNDim()); + evaluateGradient(x.data(), out.data()); + return out; +} + +std::vector RooMinimizerFcn::Hessian(std::vector const &x) const +{ + const std::size_t n = getNDim(); + std::vector out(n * n); + evaluateHessian(x, out.data()); + return out; +} + +void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer * /*context*/) { minim.SetFunction(*_multiGenFcn); - if (context->_cfg.useHessian && _funct->hasHessian()) { + if (_useHessian) { minim.SetHessianFunction( std::bind(&RooMinimizerFcn::evaluateHessian, this, std::placeholders::_1, std::placeholders::_2)); } - // The independence information for skipping vanishing second derivatives - // in numerical Hessian computations is a Minuit2-only feature, so it is - // wired up directly with the concrete minimizer type instead of going - // through the ROOT::Math::Minimizer interface. - if (auto *minuit2 = dynamic_cast(&minim)) { - minuit2->SetSecondDerivativeAlwaysVanishesFunc( - [this](unsigned int i, unsigned int j) { return secondDerivativeAlwaysVanishes(i, j); }); - } } //////////////////////////////////////////////////////////////////////////////// @@ -365,7 +377,7 @@ void RooMinimizerFcn::buildSecondDerivMask() const /// (indices in the space of all floatable parameters, matching Minuit's /// external parameter indices) is identically zero because the parameters /// share no additive term of the minimized function. -bool RooMinimizerFcn::secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const +bool RooMinimizerFcn::SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const { std::call_once(_secondDerivMaskOnce, &RooMinimizerFcn::buildSecondDerivMask, this); return !_secondDerivMask[getNDim() * i + j]; diff --git a/roofit/roofitcore/src/RooMinimizerFcn.h b/roofit/roofitcore/src/RooMinimizerFcn.h index de7052c161198..8899f8f7cf78c 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.h +++ b/roofit/roofitcore/src/RooMinimizerFcn.h @@ -49,12 +49,23 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { RooArgSet freezeDisconnectedParameters() const override; - bool secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double operator()(std::vector const &x) const override { return (*this)(x.data()); } + bool HasGradient() const override { return _useGradient; } + std::vector Gradient(std::vector const &x) const override; + bool HasHessian() const override { return _useHessian; } + std::vector Hessian(std::vector const &x) const override; + bool SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const override; + /// @} private: void buildSecondDerivMask() const; RooAbsReal *_funct = nullptr; + bool _useGradient = false; ///< Whether to provide the analytic gradient to the minimizer. + bool _useHessian = false; ///< Whether to provide the analytic Hessian to the minimizer. + /// Adapter to ROOT::Math::Minimizer, only used for minimizers other than Minuit2. std::unique_ptr _multiGenFcn; mutable std::vector _gradientOutput; mutable std::vector _hessianOutput; diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index 37b69e3acf2f3..df6995892617c 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -21,7 +21,7 @@ #include "RooMsgService.h" #include "RooMinimizer.h" -#include "Minuit2/Minuit2Minimizer.h" +#include "../RooAbsMinimizerFcn.h" #include "Minuit2/MnStrategy.h" #include @@ -193,9 +193,6 @@ void LikelihoodGradientJob::update_state() std::copy(offsets_message_begin, offsets_message_end, shared_offset_.offsets().begin()); } - // Since the gradient parallelization only support Minuit 2, we can do this cast - auto &minim = static_cast(*minimizer_->_minimizer); - // The master already knows the function value at the current point from // the line search that Minuit just completed; pre-seeding the derivator's // cache with it makes the SetupDifferentiate call below skip its full @@ -206,7 +203,7 @@ void LikelihoodGradientJob::update_state() // note: the next call must stay after the (possible) update of the offset, because it // calls the likelihood function, so the offset must be correct at this point - gradf_.SetupDifferentiate(minimizer_->getNPar(), minim.GetFCN(), minuit_internal_x_.data(), + gradf_.SetupDifferentiate(minimizer_->getNPar(), minimizer_->_fcn.get(), minuit_internal_x_.data(), minimizer_->fitter()->Config().ParamsSettings()); } } @@ -218,12 +215,9 @@ void LikelihoodGradientJob::update_state() void LikelihoodGradientJob::run_derivator(unsigned int i_component) const { - // Since the gradient parallelization only support Minuit 2, we can do this cast - auto &minim = static_cast(*minimizer_->_minimizer); - // Calculate the derivative etc for these parameters - grad_[i_component] = gradf_.FastPartialDerivative(minim.GetFCN(), minimizer_->fitter()->Config().ParamsSettings(), - i_component, grad_[i_component]); + grad_[i_component] = gradf_.FastPartialDerivative( + minimizer_->_fcn.get(), minimizer_->fitter()->Config().ParamsSettings(), i_component, grad_[i_component]); } void LikelihoodGradientJob::calculate_all() diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx index d785bfe4eabce..232312cf2ab84 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -17,57 +17,11 @@ #include "RooAbsPdf.h" #include "RooNaNPacker.h" -#include -#include - #include // std::setprecision namespace RooFit { namespace TestStatistics { -namespace { - -class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { -public: - MinuitGradFunctor(MinuitFcnGrad const &fcn, double errorLevel) : _fcn{fcn}, _up{errorLevel} {} - - double operator()(std::vector const &v) const override { return _fcn(v.data()); } - double Up() const override { return _up; } - void SetErrorDef(double val) override { _up = val; } - bool HasGradient() const override { return true; } - std::vector Gradient(std::vector const ¶ms) const override - { - std::vector grad(_fcn.getNDim()); - _fcn.Gradient(params.data(), grad.data()); - return grad; - } - // Unhide the 4-argument overload from FCNBase, which forwards to Gradient(). - // Otherwise GCC's -Woverloaded-virtual (enabled with -Werror on some CI - // targets) complains because we only override the 5-argument overload here. - using FCNBase::GradientWithPrevResult; - std::vector GradientWithPrevResult(std::vector const &v, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtV) const override - { - std::vector output(v.size()); - _fcn.GradientWithPrevResult(v.data(), output.data(), previous_grad, previous_g2, previous_gstep, fValAtV); - return output; - } - ROOT::Minuit2::GradientParameterSpace gradParameterSpace() const override - { - return _fcn.returnsInMinuit2ParameterSpace() ? ROOT::Minuit2::GradientParameterSpace::Internal - : ROOT::Minuit2::GradientParameterSpace::External; - } - - // TODO: Implement this - bool SecondDerivativeAlwaysVanishes(unsigned int /*i*/, unsigned int /*j*/) const override { return false; } - -private: - MinuitFcnGrad const &_fcn; - double _up; -}; - -} // namespace - /** \class MinuitFcnGrad * * \brief Minuit-RooMinimizer interface which synchronizes parameter data and coordinates evaluation of likelihood @@ -77,8 +31,8 @@ class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { * calculator classes to evaluate likelihood and likelihood gradient values and returns them to Minuit. The Wrapper * objects do the actual calculations. These are constructed inside the MinuitFcnGrad constructor using the RooAbsL * likelihood passed in to the constructor, usually directly from RooMinimizer, with which this class is intimately - * coupled, being a RooAbsMinimizerFcn implementation. MinuitFcnGrad inherits from ROOT::Math::IMultiGradFunction as - * well, which allows it to be used as the FCN and GRAD parameters Minuit expects. + * coupled, being a RooAbsMinimizerFcn implementation. Like all RooAbsMinimizerFcn classes, MinuitFcnGrad implements + * the ROOT::Minuit2::FCNBase interface, which allows it to be used as the FCN and GRAD parameters Minuit expects. * * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). @@ -283,10 +237,31 @@ bool MinuitFcnGrad::Synchronize(std::vector ¶m return returnee; } -void MinuitFcnGrad::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer * /*context*/) +std::vector MinuitFcnGrad::Gradient(std::vector const &x) const +{ + std::vector grad(getNDim()); + Gradient(x.data(), grad.data()); + return grad; +} + +std::vector MinuitFcnGrad::GradientWithPrevResult(std::vector const &x, double *previous_grad, + double *previous_g2, double *previous_gstep, + double fValAtX) const +{ + std::vector grad(getNDim()); + GradientWithPrevResult(x.data(), grad.data(), previous_grad, previous_g2, previous_gstep, fValAtX); + return grad; +} + +ROOT::Minuit2::GradientParameterSpace MinuitFcnGrad::gradParameterSpace() const +{ + return returnsInMinuit2ParameterSpace() ? ROOT::Minuit2::GradientParameterSpace::Internal + : ROOT::Minuit2::GradientParameterSpace::External; +} + +void MinuitFcnGrad::initMinimizer(ROOT::Math::Minimizer & /*minim*/, RooMinimizer * /*context*/) { - auto &minuit = dynamic_cast(minim); - minuit.SetFCN(getNDim(), std::make_unique(*this, minim.ErrorDef())); + throw std::logic_error("MinuitFcnGrad can only be used with Minuit2, which RooMinimizer uses directly"); } } // namespace TestStatistics diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h index f345fdf67e188..e1349bd49b1dd 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h @@ -21,7 +21,6 @@ #include "../RooAbsMinimizerFcn.h" #include -#include "Math/IFunction.h" // ROOT::Math::IMultiGradFunction class RooMinimizer; @@ -34,6 +33,8 @@ class MinuitFcnGrad : public RooAbsMinimizerFcn { std::vector ¶meters, LikelihoodMode likelihoodMode, LikelihoodGradientMode likelihoodGradientMode); + /// Not supported: this function can only be used with Minuit2, which + /// RooMinimizer drives directly via the ROOT::Minuit2::FCNBase interface. void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) override; /// Overridden from RooAbsMinimizerFcn to include gradient strategy synchronization. @@ -43,11 +44,24 @@ class MinuitFcnGrad : public RooAbsMinimizerFcn { double operator()(const double *x) const; - /// IMultiGradFunction overrides necessary for Minuit void Gradient(const double *x, double *grad) const; void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, double *previous_gstep, double fValAtX) const; + /// \name ROOT::Minuit2::FCNBase interface + /// @{ + double operator()(std::vector const &x) const override { return (*this)(x.data()); } + bool HasGradient() const override { return true; } + std::vector Gradient(std::vector const &x) const override; + // Unhide the 4-argument overload from FCNBase, which forwards to Gradient(). + // Otherwise GCC's -Woverloaded-virtual (enabled with -Werror on some CI + // targets) complains because we only override the 5-argument overload here. + using FCNBase::GradientWithPrevResult; + std::vector GradientWithPrevResult(std::vector const &x, double *previous_grad, double *previous_g2, + double *previous_gstep, double fValAtX) const override; + ROOT::Minuit2::GradientParameterSpace gradParameterSpace() const override; + /// @} + inline std::string getFunctionName() const override { return _likelihood->GetName(); } inline std::string getFunctionTitle() const override { return _likelihood->GetTitle(); } diff --git a/roofit/roofitcore/test/testRooMinimizer.cxx b/roofit/roofitcore/test/testRooMinimizer.cxx index f3f3d3cbecee1..aeaf78895987d 100644 --- a/roofit/roofitcore/test/testRooMinimizer.cxx +++ b/roofit/roofitcore/test/testRooMinimizer.cxx @@ -159,7 +159,7 @@ INSTANTIATE_TEST_SUITE_P(RooMinimizer, EvalBackendParametrizedTest, testing::Val // parameter is fixed after the RooMinimizer was constructed. The latter case // is a regression test for the translation between Minuit-internal parameter // indices (which exclude fixed parameters) and the external indices that -// RooMinimizerFcn::secondDerivativeAlwaysVanishes() is defined in. +// RooMinimizerFcn::SecondDerivativeAlwaysVanishes() is defined in. TEST(RooMinimizer, SecondDerivativeAlwaysVanishesHesse) { RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; diff --git a/roofit/xroofit/src/xRooFit.cxx b/roofit/xroofit/src/xRooFit.cxx index 6a526d549b289..f6da742d6e298 100644 --- a/roofit/xroofit/src/xRooFit.cxx +++ b/roofit/xroofit/src/xRooFit.cxx @@ -1185,7 +1185,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, }*/ // std::cout << "nIterations = " << _minimizer->fitter()->GetMinimizer()->NIterations() << std::endl; - // std::cout << "covQual before hesse = " << _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() << + // std::cout << "covQual before hesse = " << _minimizer->fitter()->Result().CovMatrixStatus() << // std::endl; sIdx = -1; if (hesseStrategy == -1) { @@ -1204,7 +1204,7 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, if (strategy == 2 && hesseStrategy == 2) { // don't repeat hesse if strategy=2 and hesseStrategy=2, and the matrix was valid - if (_minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_minimizer->fitter()->Result().CovMatrixStatus() == 3) { break; } if (sIdx >= m_hessestrategy.Length() - 1) { @@ -1260,19 +1260,19 @@ std::shared_ptr xRooFit::minimize(RooAbsReal &nll, delete _nll; throw std::runtime_error("Keyboard interrupt while hesse calculating"); } - if ((_status != 0 || _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() != 3) && status == 0 && + if ((_status != 0 || _minimizer->fitter()->Result().CovMatrixStatus() != 3) && status == 0 && printLevel >= -1) { printCerr(TString::Format("Warning: %s hesse status is %d, covQual=%d", fitName.Data(), _status, - _minimizer->fitter()->GetMinimizer()->CovMatrixStatus()) + _minimizer->fitter()->Result().CovMatrixStatus()) .Data()); } - if (_status == 0 && _minimizer->fitter()->GetMinimizer()->CovMatrixStatus() == 3) { + if (_status == 0 && _minimizer->fitter()->Result().CovMatrixStatus() == 3) { // covariance is valid! break; } else if (_status == 0) { // set the statusHistory to the cov status, since that's more informative - statusHistory.back().second = _minimizer->fitter()->GetMinimizer()->CovMatrixStatus(); + statusHistory.back().second = _minimizer->fitter()->Result().CovMatrixStatus(); } if (sIdx >= m_hessestrategy.Length() - 1) { From c1846de9a1496cc400e0ed46f7f1e5eac498e94c Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 9 Oct 2026 09:51:26 +0000 Subject: [PATCH 8/9] [Math][RF] Remove RooFit-only Minuit 2 API or move it into RooFit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Now that RooMinimizer drives Minuit2 directly, the entry points that were added to the Minuit2Minimizer wrapper only for RooFit have no users left: - Minuit2Minimizer::SetFCN(unsigned int, std::unique_ptr) and both GetFCN() overloads, which let RooFit's multiprocess gradient FCN bypass the ROOT::Math function interface. - Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(), which just forwarded to the FCNAdapter setter. RooMinimizerFcn now overrides FCNBase::SecondDerivativeAlwaysVanishes() itself. The FCNAdapter setter stays, since it is useful independently of RooFit. The NumericalDerivator class is a RooFit component that happened to live in Minuit2: nothing in Minuit2 uses it, its only user is LikelihoodGradientJob. It moves next to that user as RooFit::TestStatistics::NumericalDerivator, with the original authorship notice preserved. 🤖 Done with the help of AI --- README/ReleaseNotes/v642/index.md | 3 +- math/minuit2/CMakeLists.txt | 2 -- math/minuit2/inc/Minuit2/Minuit2Minimizer.h | 15 ---------- math/minuit2/src/Minuit2Minimizer.cxx | 15 ---------- roofit/roofitcore/CMakeLists.txt | 2 ++ .../TestStatistics/LikelihoodGradientJob.cxx | 5 ++-- .../TestStatistics/LikelihoodGradientJob.h | 9 +++--- .../TestStatistics}/NumericalDerivator.cxx | 23 ++++++++------- .../src/TestStatistics}/NumericalDerivator.h | 29 ++++++++++--------- 9 files changed, 39 insertions(+), 64 deletions(-) rename {math/minuit2/src => roofit/roofitcore/src/TestStatistics}/NumericalDerivator.cxx (97%) rename {math/minuit2/inc/Minuit2 => roofit/roofitcore/src/TestStatistics}/NumericalDerivator.h (83%) diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index a233715e4835a..09b087a24edbe 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -103,6 +103,7 @@ The `TMVA_SOFIE_GNN` tutorials have been migrated to this workflow and produce i * The header `Rstrstream.h` is deprecated and will be removed after ROOT 6.44, use instead ``. * The headers `ZipLZMA.h`, `ZipLZ4.h` and `ZipZSTD.h` are deprecated and will be removed in ROOT 6.46, use instead the public methods in the `RZip.h` interface. * Remove outdated `TSQLFile` class, was not tested for very long time, not working with only supported sqlite DB +* The `ROOT::Minuit2::Minuit2Minimizer::SetFCN(unsigned int, std::unique_ptr)` and `GetFCN()` methods are removed. They were added for RooFit, which drives Minuit2 directly now and doesn't need them anymore. To minimize an `FCNBase` object, use the Minuit2 library interface (`MnMigrad`, `MnMinimize`, etc.) directly. The `ROOT::Minuit2::NumericalDerivator` class, which was only used by RooFit, is moved to RooFit as well and is no longer part of the public Minuit2 headers. ## Build System @@ -176,7 +177,7 @@ maps) will now obtain different, mathematically consistent values. ### Skipping identically-vanishing second derivatives in numerical Hessian evaluation -Minuit2 objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function `ROOT::Minuit2::FCNBase::SecondDerivativeAlwaysVanishes()` (with corresponding setters on `ROOT::Minuit2::FCNAdapter` and `ROOT::Minuit2::Minuit2Minimizer`). +Minuit2 objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function `ROOT::Minuit2::FCNBase::SecondDerivativeAlwaysVanishes()` (with a corresponding setter on `ROOT::Minuit2::FCNAdapter`). The numerical Hessian computation in Minuit2 (`MnHesse`) skips the finite-difference evaluations for such parameter pairs, which can speed up Hesse significantly for likelihoods with many mutually independent parameters. The parameter indices in this interface always refer to the function's own full (external) parameter space, including parameters that are fixed in the minimizer, and the advertised information must hold for all parameter values; see the `FCNBase` documentation for the full contract. diff --git a/math/minuit2/CMakeLists.txt b/math/minuit2/CMakeLists.txt index aebf783361302..3e8f70b29406d 100644 --- a/math/minuit2/CMakeLists.txt +++ b/math/minuit2/CMakeLists.txt @@ -91,7 +91,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) Minuit2/ModularFunctionMinimizer.h Minuit2/NegativeG2LineSearch.h Minuit2/Numerical2PGradientCalculator.h - Minuit2/NumericalDerivator.h Minuit2/ParametricFunction.h Minuit2/ScanBuilder.h Minuit2/ScanMinimizer.h @@ -152,7 +151,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT) src/ModularFunctionMinimizer.cxx src/NegativeG2LineSearch.cxx src/Numerical2PGradientCalculator.cxx - src/NumericalDerivator.cxx src/ParametricFunction.cxx src/ScanBuilder.cxx src/SimplexBuilder.cxx diff --git a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h index c188f6038df53..7814f14ee1ed7 100644 --- a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h +++ b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h @@ -82,15 +82,6 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer { /// set the function implementing Hessian computation void SetHessianFunction(std::function, double *)> hfunc) override; - /// Set a predicate advertising which mixed second derivatives of the - /// minimized function are identically zero, so that MnHesse can skip the - /// corresponding finite-difference evaluations. The predicate must - /// fulfill the contract documented for - /// FCNBase::SecondDerivativeAlwaysVanishes(); in particular, the indices - /// refer to the full (external) parameter space. Like SetHessianFunction, - /// this must be called after SetFunction. - void SetSecondDerivativeAlwaysVanishesFunc(std::function func); - /// set free variable bool SetVariable(unsigned int ivar, const std::string &name, double val, double step) override; @@ -288,12 +279,6 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer { /// return the minimizer state (containing values, step size , etc..) const ROOT::Minuit2::MnUserParameterState &State() { return fState; } - /// To set the function directly to a Minuit 2 function. - void SetFCN(unsigned int nDim, std::unique_ptr fcn); - - const ROOT::Minuit2::FCNBase *GetFCN() const { return fMinuitFCN.get(); } - ROOT::Minuit2::FCNBase *GetFCN() { return fMinuitFCN.get(); } - protected: // protected function for accessing the internal Minuit2 object. Needed for derived classes diff --git a/math/minuit2/src/Minuit2Minimizer.cxx b/math/minuit2/src/Minuit2Minimizer.cxx index 54b87884ae7b1..94e401325e035 100644 --- a/math/minuit2/src/Minuit2Minimizer.cxx +++ b/math/minuit2/src/Minuit2Minimizer.cxx @@ -416,15 +416,6 @@ void Minuit2Minimizer::SetHessianFunction(std::functionSetHessianFunction(hfunc); } -void Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(std::function func) -{ - // not supported for Fumili, whose FCN is not an FCNAdapter - if (fUseFumili) return; - auto fcn = static_cast(fMinuitFCN.get()); - if (!fcn) return; - fcn->SetSecondDerivativeAlwaysVanishesFunc(std::move(func)); -} - namespace { ROOT::Minuit2::MnStrategy customizedStrategy(unsigned int strategyLevel, ROOT::Math::MinimizerOptions const &options) @@ -1372,12 +1363,6 @@ void Minuit2Minimizer::SetStorageLevel(int level) fMinimizer->Builder().SetStorageLevel(level); } -void Minuit2Minimizer::SetFCN(unsigned int nDim, std::unique_ptr fcn) -{ - fDim = nDim; - fMinuitFCN = std::move(fcn); -} - } // end namespace Minuit2 } // end namespace ROOT diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 02527b41db31b..3cef3d2e82b80 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -14,12 +14,14 @@ if(NOT WIN32) src/TestStatistics/LikelihoodGradientJob.cxx src/TestStatistics/LikelihoodJob.cxx src/TestStatistics/MinuitFcnGrad.cxx + src/TestStatistics/NumericalDerivator.cxx ) set(RooFitMPTestStatisticsHeaders src/TestStatistics/LikelihoodGradientJob.h src/TestStatistics/LikelihoodSerial.h src/TestStatistics/MinuitFcnGrad.h src/TestStatistics/LikelihoodJob.h + src/TestStatistics/NumericalDerivator.h ) list(APPEND EXTRA_LIBRARIES RooFitMultiProcess) #FIXME: The ProcessTimer.h exposes json in its interface: diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index df6995892617c..b66f424853e93 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -169,9 +169,8 @@ void LikelihoodGradientJob::update_state() auto gradient_message = get_manager()->messenger().receive_from_master_on_worker(&more); assert(more); - auto gradient_message_begin = gradient_message.data(); - auto gradient_message_end = - gradient_message_begin + gradient_message.size() / sizeof(ROOT::Minuit2::DerivatorElement); + auto gradient_message_begin = gradient_message.data(); + auto gradient_message_end = gradient_message_begin + gradient_message.size() / sizeof(DerivatorElement); std::copy(gradient_message_begin, gradient_message_end, grad_.begin()); auto minuit_internal_x_message = diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h index da94b01efb800..5951d64b8eb8d 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h @@ -17,9 +17,10 @@ #include "RooFit/TestStatistics/LikelihoodGradientWrapper.h" #include "Math/MinimizerOptions.h" -#include "Minuit2/NumericalDerivator.h" #include "Minuit2/MnMatrix.h" +#include "NumericalDerivator.h" + #include #include @@ -63,7 +64,7 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien struct task_result_t { std::size_t job_id; std::size_t task_id; - ROOT::Minuit2::DerivatorElement grad; + DerivatorElement grad; }; void send_back_task_result_from_worker(std::size_t task) override; bool receive_task_result_on_master(const RooFit::MultiProcess::Message &message) override; @@ -75,8 +76,8 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien // members // mutables below are because ROOT::Math::IMultiGradFunction::DoDerivative is const - mutable std::vector grad_; - mutable ROOT::Minuit2::NumericalDerivator gradf_; + mutable std::vector grad_; + mutable NumericalDerivator gradf_; std::size_t N_tasks_ = 0; std::size_t N_tasks_at_workers_ = 0; diff --git a/math/minuit2/src/NumericalDerivator.cxx b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx similarity index 97% rename from math/minuit2/src/NumericalDerivator.cxx rename to roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx index d028ed95d7c84..5151cf44116bd 100644 --- a/math/minuit2/src/NumericalDerivator.cxx +++ b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.cxx @@ -1,4 +1,3 @@ -// @(#)root/mathcore:$Id$ // Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 // /********************************************************************** @@ -23,18 +22,20 @@ * Modified things (w.r.t. original) are indicated by MODIFIED. */ -#include "Minuit2/NumericalDerivator.h" -#include -#include -#include +#include "NumericalDerivator.h" + +#include #include + +#include #include -#include "Fit/ParameterSettings.h" +#include +#include -#include // needed here because in Fitter is only a forward declaration +namespace RooFit { +namespace TestStatistics { -namespace ROOT { -namespace Minuit2 { +using ROOT::Minuit2::FCNBase; NumericalDerivator::NumericalDerivator(bool always_exactly_mimic_minuit2) : fAlwaysExactlyMimicMinuit2(always_exactly_mimic_minuit2) @@ -286,5 +287,5 @@ std::ostream &operator<<(std::ostream &out, const DerivatorElement &value) << ", step_size: " << value.step_size << ")"; } -} // namespace Minuit2 -} // namespace ROOT +} // namespace TestStatistics +} // namespace RooFit diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.h similarity index 83% rename from math/minuit2/inc/Minuit2/NumericalDerivator.h rename to roofit/roofitcore/src/TestStatistics/NumericalDerivator.h index f8e42fc35fd7a..ad2c83001781a 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/roofit/roofitcore/src/TestStatistics/NumericalDerivator.h @@ -1,4 +1,3 @@ -// @(#)root/mathcore:$Id$ // Authors: L. Moneta, J.T. Offermann, E.G.P. Bos 2013-2018 // /********************************************************************** @@ -13,10 +12,13 @@ * Authors: L. Moneta, J. T. Offermann * Modified version created on: Sep 27, 2017 * Author: E. G. P. Bos + * + * Lived in Minuit2 as ROOT::Minuit2::NumericalDerivator until it was moved + * to RooFit, its only user. */ -#ifndef ROOT_Minuit2_NumericalDerivator -#define ROOT_Minuit2_NumericalDerivator +#ifndef ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator +#define ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator #include "Fit/ParameterSettings.h" #include "Minuit2/MnParameterTransformation.h" @@ -27,8 +29,8 @@ #include -namespace ROOT { -namespace Minuit2 { +namespace RooFit { +namespace TestStatistics { // Holds all necessary derivatives and associated numbers (per parameter) used in the NumericalDerivator class. struct DerivatorElement { @@ -44,7 +46,7 @@ class NumericalDerivator { NumericalDerivator(double step_tolerance, double grad_tolerance, unsigned int ncycles, double error_level, bool always_exactly_mimic_minuit2 = true); - void SetupDifferentiate(unsigned int nDim, const FCNBase *function, const double *cx, + void SetupDifferentiate(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *cx, std::span parameters); /// Pre-seed the cache of the function value at the central point, so that a @@ -58,17 +60,18 @@ class NumericalDerivator { fVal = fval; } - std::vector Differentiate(unsigned int nDim, const FCNBase *function, const double *x, + std::vector Differentiate(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, + const double *x, std::span parameters, std::span previous_gradient); - DerivatorElement PartialDerivative(unsigned int nDim, const FCNBase *function, const double *x, + DerivatorElement PartialDerivative(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *x, std::span parameters, unsigned int i_component, DerivatorElement previous); - DerivatorElement FastPartialDerivative(const FCNBase *function, + DerivatorElement FastPartialDerivative(const ROOT::Minuit2::FCNBase *function, std::span parameters, unsigned int i_component, const DerivatorElement &previous); - DerivatorElement operator()(unsigned int nDim, const FCNBase *function, const double *x, + DerivatorElement operator()(unsigned int nDim, const ROOT::Minuit2::FCNBase *function, const double *x, std::span parameters, unsigned int i_component, const DerivatorElement &previous); @@ -114,7 +117,7 @@ class NumericalDerivator { std::ostream &operator<<(std::ostream &out, const DerivatorElement &value); -} // namespace Minuit2 -} // namespace ROOT +} // namespace TestStatistics +} // namespace RooFit -#endif // ROOT_Minuit2_NumericalDerivator +#endif // ROOT_ROOFIT_TESTSTATISTICS_NumericalDerivator From 09938e17c9e401c2acfda198c4fbc1803ae93078 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 9 Oct 2026 10:06:26 +0000 Subject: [PATCH 9/9] [RF] Mechanical cleanup after driving Minuit2 directly from RooMinimizer MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A collection of small simplifications that fall out of RooMinimizer using Minuit2 directly and the FCN classes implementing ROOT::Minuit2::FCNBase: - RooMinimizerFcn::initMinimizer() no longer passes the Hessian through ROOT::Math::Minimizer::SetHessianFunction(). Only the Minuit2 wrapper implements that method, and Minuit2 doesn't go through this path anymore, so it was a no-op. The Hessian is now computed directly in the FCNBase::Hessian() override. - initMinimizer() is a plain member of RooMinimizerFcn instead of a virtual on RooAbsMinimizerFcn, since the multiprocess MinuitFcnGrad could only throw there. RooMinimizer checks the function type when it creates a ROOT::Math::Minimizer. - The raw-pointer evaluation and gradient methods of MinuitFcnGrad were only wrappers called by the FCNBase overrides; the implementation moves into the overrides. - The RooMinimizer::FitResult constructor from a FitConfig and the fallback that used it in calculateHessErrors() were unreachable: fitFCN() rejects functions without parameters before a result is created, and hesse() refuses to run without a result. - The FCN classes and LikelihoodGradientJob access the RooMinimizer members directly instead of going through the fitter() compatibility shim, which allocated a FitterInterface on every call. 🤖 Done with the help of AI --- roofit/roofitcore/inc/RooMinimizer.h | 1 - roofit/roofitcore/src/RooAbsMinimizerFcn.cxx | 2 +- roofit/roofitcore/src/RooAbsMinimizerFcn.h | 8 +--- roofit/roofitcore/src/RooMinimizer.cxx | 43 ++++------------- roofit/roofitcore/src/RooMinimizerFcn.cxx | 19 ++------ roofit/roofitcore/src/RooMinimizerFcn.h | 8 +++- .../TestStatistics/LikelihoodGradientJob.cxx | 6 +-- .../src/TestStatistics/MinuitFcnGrad.cxx | 48 +++++++------------ .../src/TestStatistics/MinuitFcnGrad.h | 12 +---- .../testLikelihoodGradientJob.cxx | 6 +-- 10 files changed, 45 insertions(+), 108 deletions(-) diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 2238974c3bbfe..15ef5b00af88e 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -49,7 +49,6 @@ class RooMinimizer : public TObject { struct FitResult { FitResult() = default; - FitResult(const ROOT::Fit::FitConfig &fconfig); double error(unsigned int i) const { return (i < fErrors.size()) ? fErrors[i] : 0; } double lowerError(unsigned int i) const; diff --git a/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx b/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx index 1e5ed00c84593..bc35a91f3af20 100644 --- a/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx @@ -154,7 +154,7 @@ bool RooAbsMinimizerFcn::Synchronize(std::vector & /// Transfer MINUIT fit results back into RooFit objects. void RooAbsMinimizerFcn::BackProp() { - auto const &results = _context->fitter()->Result(); + auto const &results = *_context->_result; for (std::size_t index = 0; index < getNDim(); index++) { diff --git a/roofit/roofitcore/src/RooAbsMinimizerFcn.h b/roofit/roofitcore/src/RooAbsMinimizerFcn.h index 3c2fe7e8730b0..cb91ad7196e0c 100644 --- a/roofit/roofitcore/src/RooAbsMinimizerFcn.h +++ b/roofit/roofitcore/src/RooAbsMinimizerFcn.h @@ -22,7 +22,7 @@ #include "RooMinimizer.h" #include "RooRealVar.h" -#include +#include #include #include @@ -41,11 +41,6 @@ class RooAbsMinimizerFcn : public ROOT::Minuit2::FCNBase { RooAbsMinimizerFcn(RooArgList paramList, RooMinimizer *context); ~RooAbsMinimizerFcn() override = default; - /// Set this function on a ROOT::Math::Minimizer. Only used for minimizers - /// other than Minuit2, which RooMinimizer uses directly via the FCNBase - /// interface that this class implements. - virtual void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) = 0; - /// \name ROOT::Minuit2::FCNBase interface /// @{ double Up() const override { return _errorDef; } @@ -94,6 +89,7 @@ class RooAbsMinimizerFcn : public ROOT::Minuit2::FCNBase { virtual void setOffsetting(bool flag) = 0; RooMinimizer::Config const &cfg() const { return _context->_cfg; } + ROOT::Math::MinimizerOptions const &minimizerOptions() const { return _context->_config.MinimizerOptions(); } inline RooRealVar &floatableParam(std::size_t i) const { diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 9f32ff6ce674c..0df6fff8999e1 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -1232,11 +1232,6 @@ bool RooMinimizer::calculateHessErrors() if (!ret) coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error when calculating Hessian" << std::endl; - // update minimizer results with what comes out from Hesse - // in case is empty - create from a FitConfig - if (_result->fParams.empty()) - _result = std::make_unique(_config); - // re-give a minimizer instance in case it has been changed ret |= update(ret); @@ -1320,8 +1315,16 @@ void RooMinimizer::initMinimizer() initMinuit2(); return; } + // Only RooMinimizerFcn supports minimizers other than Minuit2: the + // multiprocess MinuitFcnGrad is tied to Minuit2, which the constructor + // enforces. + auto *fcn = dynamic_cast(_fcn.get()); + if (!fcn) { + throw std::logic_error("RooMinimizer: minimizer type " + _cfg.minimizerType + + " is only supported with the default (non-parallel) likelihood function"); + } _minimizer = std::unique_ptr(_config.CreateMinimizer()); - _fcn->initMinimizer(*_minimizer, this); + fcn->initMinimizer(*_minimizer); _minimizer->SetVariables(_config.ParamsSettings().begin(), _config.ParamsSettings().end()); if (_cfg.setInitialCovariance) { @@ -1376,34 +1379,6 @@ void RooMinimizer::updateFitConfig() } } -RooMinimizer::FitResult::FitResult(const ROOT::Fit::FitConfig &fconfig) - : fStatus(-99), // use this special convention to flag it when printing result - fCovStatus(0), - fParams(fconfig.NPar()), - fErrors(fconfig.NPar()) -{ - // create a Fit result from a fit config (i.e. with initial parameter values - // and errors equal to step values - // The model function is NULL in this case - - // set minimizer type and algorithm - fMinimType = fconfig.MinimizerType(); - // append algorithm name for minimizer that support it - if ((fMinimType.find("Fumili") == std::string::npos) && (fMinimType.find("GSLMultiFit") == std::string::npos)) { - if (!fconfig.MinimizerAlgoType().empty()) - fMinimType += " / " + fconfig.MinimizerAlgoType(); - } - - // get parameter values and errors (step sizes) - for (unsigned int i = 0; i < fconfig.NPar(); ++i) { - const ROOT::Fit::ParameterSettings &par = fconfig.ParSettings(i); - fParams[i] = par.Value(); - fErrors[i] = par.StepSize(); - if (par.IsFixed()) - fFixedParams[i] = true; - } -} - void RooMinimizer::fillResult(bool isValid) { if (useMinuit2Directly()) { diff --git a/roofit/roofitcore/src/RooMinimizerFcn.cxx b/roofit/roofitcore/src/RooMinimizerFcn.cxx index efe06b57e8184..3f482f65b515f 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooMinimizerFcn.cxx @@ -272,7 +272,7 @@ RooArgSet RooMinimizerFcn::freezeDisconnectedParameters() const return changedSet; } -bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) const +std::vector RooMinimizerFcn::Hessian(std::vector const &x) const { // Set the parameter values for this iteration for (unsigned index = 0; index < getNDim(); index++) { @@ -285,6 +285,7 @@ bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) co std::size_t m = _allParamsInit.size(); std::size_t n = getNDim(); + std::vector out(n * n); std::size_t iAll = 0; std::size_t iFloating = 0; for (RooAbsArg *param_i : _allParamsInit) { @@ -313,7 +314,7 @@ bool RooMinimizerFcn::evaluateHessian(std::span x, double *out) co std::cout << std::endl; } } - return true; + return out; } std::vector RooMinimizerFcn::Gradient(std::vector const &x) const @@ -323,21 +324,9 @@ std::vector RooMinimizerFcn::Gradient(std::vector const &x) cons return out; } -std::vector RooMinimizerFcn::Hessian(std::vector const &x) const -{ - const std::size_t n = getNDim(); - std::vector out(n * n); - evaluateHessian(x, out.data()); - return out; -} - -void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer * /*context*/) +void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim) const { minim.SetFunction(*_multiGenFcn); - if (_useHessian) { - minim.SetHessianFunction( - std::bind(&RooMinimizerFcn::evaluateHessian, this, std::placeholders::_1, std::placeholders::_2)); - } } //////////////////////////////////////////////////////////////////////////////// diff --git a/roofit/roofitcore/src/RooMinimizerFcn.h b/roofit/roofitcore/src/RooMinimizerFcn.h index 8899f8f7cf78c..760a9015d9abb 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.h +++ b/roofit/roofitcore/src/RooMinimizerFcn.h @@ -22,6 +22,8 @@ #include "RooAbsReal.h" #include "RooArgList.h" +#include + #include #include #include @@ -36,7 +38,10 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { public: RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context); - void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) override; + /// Set this function on a ROOT::Math::Minimizer. Only used for minimizers + /// other than Minuit2, which RooMinimizer drives directly via the + /// ROOT::Minuit2::FCNBase interface. + void initMinimizer(ROOT::Math::Minimizer &minim) const; std::string getFunctionName() const override; std::string getFunctionTitle() const override; @@ -45,7 +50,6 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { double operator()(const double *x) const; void evaluateGradient(const double *x, double *out) const; - bool evaluateHessian(std::span x, double *out) const; RooArgSet freezeDisconnectedParameters() const override; diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index b66f424853e93..242d37e72482d 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -203,7 +203,7 @@ void LikelihoodGradientJob::update_state() // note: the next call must stay after the (possible) update of the offset, because it // calls the likelihood function, so the offset must be correct at this point gradf_.SetupDifferentiate(minimizer_->getNPar(), minimizer_->_fcn.get(), minuit_internal_x_.data(), - minimizer_->fitter()->Config().ParamsSettings()); + minimizer_->_config.ParamsSettings()); } } @@ -215,8 +215,8 @@ void LikelihoodGradientJob::update_state() void LikelihoodGradientJob::run_derivator(unsigned int i_component) const { // Calculate the derivative etc for these parameters - grad_[i_component] = gradf_.FastPartialDerivative( - minimizer_->_fcn.get(), minimizer_->fitter()->Config().ParamsSettings(), i_component, grad_[i_component]); + grad_[i_component] = gradf_.FastPartialDerivative(minimizer_->_fcn.get(), minimizer_->_config.ParamsSettings(), + i_component, grad_[i_component]); } void LikelihoodGradientJob::calculate_all() diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx index 232312cf2ab84..1fdfc890da782 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -95,9 +95,9 @@ void MinuitFcnGrad::syncOffsets() const } } -double MinuitFcnGrad::operator()(const double *x) const +double MinuitFcnGrad::operator()(std::vector const &x) const { - syncParameterValuesFromMinuitCalls(x, false); + syncParameterValuesFromMinuitCalls(x.data(), false); syncOffsets(); @@ -207,23 +207,28 @@ bool MinuitFcnGrad::syncParameterValuesFromMinuitCalls(const double *x, bool min return aParamWasUpdated; } -void MinuitFcnGrad::Gradient(const double *x, double *grad) const +std::vector MinuitFcnGrad::Gradient(std::vector const &x) const { + std::vector grad(getNDim()); _calculatingGradient = true; - syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); + syncParameterValuesFromMinuitCalls(x.data(), returnsInMinuit2ParameterSpace()); syncOffsets(); - _gradient->fillGradient(grad); + _gradient->fillGradient(grad.data()); _calculatingGradient = false; + return grad; } -void MinuitFcnGrad::GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtX) const +std::vector MinuitFcnGrad::GradientWithPrevResult(std::vector const &x, double *previous_grad, + double *previous_g2, double *previous_gstep, + double fValAtX) const { + std::vector grad(getNDim()); _calculatingGradient = true; - syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); + syncParameterValuesFromMinuitCalls(x.data(), returnsInMinuit2ParameterSpace()); syncOffsets(); - _gradient->fillGradientWithPrevResult(grad, previous_grad, previous_g2, previous_gstep, fValAtX); + _gradient->fillGradientWithPrevResult(grad.data(), previous_grad, previous_g2, previous_gstep, fValAtX); _calculatingGradient = false; + return grad; } bool MinuitFcnGrad::Synchronize(std::vector ¶meters) @@ -232,37 +237,16 @@ bool MinuitFcnGrad::Synchronize(std::vector ¶m applyToLikelihood([&](auto &l) { l.synchronizeParameterSettings(parameters); }); _gradient->synchronizeParameterSettings(parameters); - applyToLikelihood([&](auto &l) { l.synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); }); - _gradient->synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); + applyToLikelihood([&](auto &l) { l.synchronizeWithMinimizer(minimizerOptions()); }); + _gradient->synchronizeWithMinimizer(minimizerOptions()); return returnee; } -std::vector MinuitFcnGrad::Gradient(std::vector const &x) const -{ - std::vector grad(getNDim()); - Gradient(x.data(), grad.data()); - return grad; -} - -std::vector MinuitFcnGrad::GradientWithPrevResult(std::vector const &x, double *previous_grad, - double *previous_g2, double *previous_gstep, - double fValAtX) const -{ - std::vector grad(getNDim()); - GradientWithPrevResult(x.data(), grad.data(), previous_grad, previous_g2, previous_gstep, fValAtX); - return grad; -} - ROOT::Minuit2::GradientParameterSpace MinuitFcnGrad::gradParameterSpace() const { return returnsInMinuit2ParameterSpace() ? ROOT::Minuit2::GradientParameterSpace::Internal : ROOT::Minuit2::GradientParameterSpace::External; } -void MinuitFcnGrad::initMinimizer(ROOT::Math::Minimizer & /*minim*/, RooMinimizer * /*context*/) -{ - throw std::logic_error("MinuitFcnGrad can only be used with Minuit2, which RooMinimizer uses directly"); -} - } // namespace TestStatistics } // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h index e1349bd49b1dd..7e62676456ce8 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h @@ -33,24 +33,14 @@ class MinuitFcnGrad : public RooAbsMinimizerFcn { std::vector ¶meters, LikelihoodMode likelihoodMode, LikelihoodGradientMode likelihoodGradientMode); - /// Not supported: this function can only be used with Minuit2, which - /// RooMinimizer drives directly via the ROOT::Minuit2::FCNBase interface. - void initMinimizer(ROOT::Math::Minimizer &, RooMinimizer *context) override; - /// Overridden from RooAbsMinimizerFcn to include gradient strategy synchronization. bool Synchronize(std::vector ¶meter_settings) override; bool returnsInMinuit2ParameterSpace() const { return _gradient->usesMinuitInternalValues(); } - double operator()(const double *x) const; - - void Gradient(const double *x, double *grad) const; - void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep, double fValAtX) const; - /// \name ROOT::Minuit2::FCNBase interface /// @{ - double operator()(std::vector const &x) const override { return (*this)(x.data()); } + double operator()(std::vector const &x) const override; bool HasGradient() const override { return true; } std::vector Gradient(std::vector const &x) const override; // Unhide the 4-argument overload from FCNBase, which forwards to Gradient(). diff --git a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx index db0641eba7977..5f1f44785ae3f 100644 --- a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx +++ b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx @@ -968,10 +968,10 @@ TEST(MinuitFcnGrad, DISABLED_CompareToRooMinimizerFcn) RFTS::LikelihoodGradientMode::multiprocess); RooMinimizerFcn vanilla_fcn(nll_vanilla.get(), &m_vanilla); - EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn).data()), modularL_fcn(getParamVals(modularL_fcn).data())); + EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn)), modularL_fcn(getParamVals(modularL_fcn))); // let's also check with absolutely certain same parameter values, both of them - EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn).data()), modularL_fcn(getParamVals(vanilla_fcn).data())); - EXPECT_EQ(vanilla_fcn(getParamVals(modularL_fcn).data()), modularL_fcn(getParamVals(modularL_fcn).data())); + EXPECT_EQ(vanilla_fcn(getParamVals(vanilla_fcn)), modularL_fcn(getParamVals(vanilla_fcn))); + EXPECT_EQ(vanilla_fcn(getParamVals(modularL_fcn)), modularL_fcn(getParamVals(modularL_fcn))); // reset static variables to automatic RFMP::Config::LikelihoodJob::defaultNEventTasks = RFMP::Config::LikelihoodJob::automaticNEventTasks;