[RF] Test RooFit Hessians with Clad

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}

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();

math/mathcore/inc/Math/Functor.h

@@ -23,6 +23,7 @@
 #include <algorithm>
 #include <memory>
 #include <functional>
+#include <stdexcept>
 #include <type_traits>
 #include <vector>
 
@@ -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<double(double const *)> const&f, unsigned int dim,
-               std::function<void(double const *, double *)> const& g)
-      : fDim{dim}, fFunc{f}, fGradFunc{g}
+   GradFunctor(std::function<double(double const *)> const &f, unsigned int dim,
+               std::function<void(double const *, double *)> const &g,
+               std::function<void(double const *, double *)> 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<double(const double *)> fFunc;
    std::function<double(double const *, unsigned int)> fDerivFunc;
    std::function<void(const double *, double*)> fGradFunc;
+   std::function<void(const double *, double *)> 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 <class PtrObj, typename MemFn, typename GradMemFn>
-   GradFunctor1D(const PtrObj& p, MemFn memFn, GradMemFn gradFn)
+   template <class PtrObj, typename MemFn, typename GradMemFn,
+             std::enable_if_t<std::is_member_pointer<MemFn>::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<double(double)> const& f, std::function<double(double)> 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<double(double)> const &f, std::function<double(double)> const &g,
+                 std::function<double(double)> 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<double(double)> fFunc;
    std::function<double(double)> fDerivFunc;
+   std::function<double(double)> fSecondDerivFunc;
 };
 
 

math/mathcore/inc/Math/IFunction.h

@@ -33,6 +33,7 @@
 
 #include "Math/IFunctionfwd.h"
 
+#include <stdexcept>
 
 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())");
+         }
       };
 
 

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 <plugins/include/clad/Differentiator/BuiltinDerivatives.h>\n"
+                 "namespace clad::custom_derivatives {\n\n"
+                 "clad::ValueAndPushforward<double, double> " +
+                 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<double, "
+                 "double> 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());
    }
 
@@ -590,7 +651,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()));
 }
 

roofit/histfactory/test/testHistFactory.cxx

@@ -12,17 +12,19 @@
 #include <RooFitHS3/JSONIO.h>
 #include <RooFitHS3/RooJSONFactoryWSTool.h>
 
-#include <RooFit/Detail/NormalizationHelpers.h>
-#include <RooDataHist.h>
-#include <RooWorkspace.h>
 #include <RooArgSet.h>
-#include <RooSimultaneous.h>
-#include <RooRealSumPdf.h>
-#include <RooRealVar.h>
-#include <RooHelpers.h>
+#include <RooDataHist.h>
+#include <RooEvaluatorWrapper.h>
+#include <RooFit/Detail/NormalizationHelpers.h>
+#include <RooFit/Evaluator.h>
 #include <RooFitResult.h>
+#include <RooHelpers.h>
+#include <RooMinimizer.h>
 #include <RooPlot.h>
-#include <RooFit/Evaluator.h>
+#include <RooRealSumPdf.h>
+#include <RooRealVar.h>
+#include <RooSimultaneous.h>
+#include <RooWorkspace.h>
 
 #include <TROOT.h>
 #include <TFile.h>
@@ -610,8 +612,21 @@ TEST_P(HFFixtureFit, Fit)
    }
 
    using namespace RooFit;
-   std::unique_ptr<RooFitResult> fitResult{simPdf->fitTo(
-      *data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()), Save(), PrintLevel(verbose ? 1 : -1))};
+   std::unique_ptr<RooAbsReal> 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<RooFit::Experimental::RooEvaluatorWrapper &>(*nll).generateGradient();
+      static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*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<RooFitResult> fitResult{minim.save()};
    ASSERT_NE(fitResult, nullptr);
    if (verbose)
       fitResult->Print("v");