From c9f57cf85b3761ea02336563e48da7e230af4c09 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 11 Aug 2026 12:16:55 +0000 Subject: [PATCH 1/8] [RF] Fix the codegen backend for the standard-parametrized RooLognormal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The code generated for a RooLognormal with useStandardParametrization() called RooFit::Detail::MathFuncs::logNormalEvaluateStandard, which does not exist -- the function is called logNormalStandard. Any codegen or AD fit of such a pdf failed to compile with "no member named 'logNormalEvaluateStandard'". This went unnoticed because the LognormalStandard case in testRooFuncWrapper builds its pdf with Lognormal::model(x[...], mu[...], k[...], true) and the factory quietly dropped that last argument, so the test was really a duplicate of the Lognormal one. RooFactoryWSTool::asINT(), which is also the conversion used for bool constructor arguments, is atoi(), and atoi("true") is zero. Teach it about the spelled-out literals, which affects every factory string that writes a bool that way, and makes the existing test exercise what it says it does. Verified that the generated code now agrees exactly with the reference backend, for the nominal likelihood value and over a scan of the shape parameters. 🤖 Done with the help of AI --- roofit/codegen/src/CodegenImpl.cxx | 2 +- roofit/roofitcore/src/RooFactoryWSTool.cxx | 9 ++++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index 750f6912bc380..15598eb27b96b 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -590,7 +590,7 @@ void codegenImpl(RooLandau &arg, CodegenContext &ctx) void codegenImpl(RooLognormal &arg, CodegenContext &ctx) { - std::string funcName = arg.useStandardParametrization() ? "logNormalEvaluateStandard" : "logNormal"; + std::string funcName = arg.useStandardParametrization() ? "logNormalStandard" : "logNormal"; ctx.addResult(&arg, ctx.buildCall(mathFunc(funcName), arg.getX(), arg.getShapeK(), arg.getMedian())); } diff --git a/roofit/roofitcore/src/RooFactoryWSTool.cxx b/roofit/roofitcore/src/RooFactoryWSTool.cxx index adcd24acf0889..ef36fb8b71188 100644 --- a/roofit/roofitcore/src/RooFactoryWSTool.cxx +++ b/roofit/roofitcore/src/RooFactoryWSTool.cxx @@ -1873,7 +1873,14 @@ const char* RooFactoryWSTool::asSTRING(const char* arg) Int_t RooFactoryWSTool::asINT(const char* arg) { - return atoi(arg) ; + // This is also the conversion used for bool constructor arguments, where + // atoi() would silently turn the spelled-out literals into zero. + const std::string s{arg}; + if (s == "true") + return 1; + if (s == "false") + return 0; + return atoi(arg); } From ee9dfee6c7ee83d2919263df84f8c5fc7f5b0697 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 16 Mar 2026 18:53:11 +0100 Subject: [PATCH 2/8] [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 9a5250e8b25b7..88b390703da3f 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -129,14 +129,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()}; } @@ -164,6 +166,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 54ac2d00479702d7d09372749d660e154f38c687 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:45:48 +0000 Subject: [PATCH 3/8] [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 53645ed5f253ba0bd7f6723113736c9e9dad929a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:46:04 +0000 Subject: [PATCH 4/8] [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 15598eb27b96b..574fec2480a0c 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -483,6 +483,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 98720b3aa79b141fcbe7f40fb5c0853262f2c0e7 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 13:46:21 +0000 Subject: [PATCH 5/8] [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 | 194 ++++++++++++++---- 1 file changed, 152 insertions(+), 42 deletions(-) diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index 88b390703da3f..91095a201fa4f 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) + double fitResultTolerance, bool randomizeParameters, double hesseTolerance = -1., + bool enableHessian = true) : _name{name}, _setupWorkspace{setupWorkspace}, _createNLL{createNLL}, _fitResultTolerance{fitResultTolerance}, - _randomizeParameters{randomizeParameters} + _hesseTolerance{hesseTolerance}, + _randomizeParameters{randomizeParameters}, + _enableHessian{enableHessian} { } @@ -110,7 +115,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; }; class FactoryTest : public testing::TestWithParam { @@ -129,11 +142,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); @@ -167,10 +180,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); @@ -207,7 +222,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 @@ -219,8 +234,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. @@ -228,7 +243,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)"); @@ -239,7 +257,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", @@ -274,7 +293,8 @@ FactoryTestParams param4{"ConstraintSum", pdf.createNLL(data, ExternalConstraints(*ws.pdf("fconstext")), backend)}; }, 1e-4, - /*randomizeParameters=*/true}; + /*randomizeParameters=*/true, + /*hesseTolerance=*/1e-3}; namespace { @@ -346,12 +366,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) { @@ -393,12 +415,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) { @@ -409,7 +433,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) { @@ -421,7 +446,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) { @@ -432,9 +458,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(); @@ -449,7 +476,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) { @@ -508,10 +536,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; } @@ -535,9 +567,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; } @@ -549,10 +607,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; } @@ -561,12 +623,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}); @@ -577,15 +644,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) + double fitResultTolerance, bool randomizeParameters = true, + double hesseTolerance = -1., bool enableHessian = true, std::size_t nEvents = 0) { return {name, [=](RooWorkspace &ws) { @@ -593,12 +672,23 @@ FactoryTestParams makeTestParams(const char *name, std::vector cons ws.factory(expr.c_str()); } ws.defineSet("observables", "x"); + // Provide a custom dataset if the default number of events that + // the test generates otherwise is not right for this model. + if (nEvents > 0) { + std::unique_ptr data{ws.pdf("model")->generate(*ws.set("observables"), nEvents)}; + std::unique_ptr binned{data->binnedClone()}; + binned->SetName("data"); + ws.import(*binned); + } }, [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) { using namespace RooFit; return std::unique_ptr{pdf.createNLL(data, backend)}; }, - fitResultTolerance, randomizeParameters}; + fitResultTolerance, + randomizeParameters, + hesseTolerance, + enableHessian}; } auto testValues = testing::Values( @@ -606,12 +696,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, -10, 10]})", + {"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", @@ -625,22 +718,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), - makeTestParams("RooLandau3", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[10., 0.01, 50.])"}, 7e-3, false), - 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), + 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, + /*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 657692e977a3d9c61732819bf2ced1f4850ab6b6 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 15:04:41 +0000 Subject: [PATCH 6/8] [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 342eb2be75858aaac075bc06e92fb960ffa08462 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 12 Aug 2026 21:12:03 +0000 Subject: [PATCH 7/8] [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 14475ea61d73223cd5c8b06f2355a6d9c794ac35 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 11 Aug 2026 12:01:17 +0200 Subject: [PATCH 8/8] Use a Clad branch that should have all the Hessian fixes --- interpreter/cling/tools/plugins/clad/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/interpreter/cling/tools/plugins/clad/CMakeLists.txt b/interpreter/cling/tools/plugins/clad/CMakeLists.txt index ab93259e12fa1..0298aad36a1f2 100644 --- a/interpreter/cling/tools/plugins/clad/CMakeLists.txt +++ b/interpreter/cling/tools/plugins/clad/CMakeLists.txt @@ -71,7 +71,7 @@ if (DEFINED CLAD_SOURCE_DIR) list(APPEND _clad_extra_settings SOURCE_DIR ${CLAD_SOURCE_DIR}) else() list(APPEND _clad_extra_settings GIT_REPOSITORY https://github.com/vgvassilev/clad.git) - list(APPEND _clad_extra_settings GIT_TAG v2.4) + list(APPEND _clad_extra_settings GIT_TAG master) endif() # list(APPEND _clad_patches_list "patch1.patch" "patch2.patch") @@ -86,7 +86,7 @@ set(_clad_patch_command ExternalProject_Add( clad UPDATE_COMMAND "" - PATCH_COMMAND ${_clad_patch_command} + # PATCH_COMMAND ${_clad_patch_command} CMAKE_ARGS -G ${CMAKE_GENERATOR} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}