[RF] Test RooFit Hessians with Clad

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 `<sstream>`.
 * 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<FCNBase>)` 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.
 

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())");
+         }
       };
 
 

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

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<bool(std::span<const double>, 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<bool(unsigned int, unsigned int)> 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<ROOT::Minuit2::FCNBase> 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
 

math/minuit2/src/Minuit2Minimizer.cxx

@@ -416,15 +416,6 @@ void Minuit2Minimizer::SetHessianFunction(std::function<bool(std::span<const dou
    fcn->SetHessianFunction(hfunc);
 }
 
-void Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(std::function<bool(unsigned int, unsigned int)> func)
-{
-   // not supported for Fumili, whose FCN is not an FCNAdapter
-   if (fUseFumili) return;
-   auto fcn = static_cast<ROOT::Minuit2::FCNAdapter *>(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<ROOT::Minuit2::FCNBase> fcn)
-{
-   fDim = nDim;
-   fMinuitFCN = std::move(fcn);
-}
-
 } // end namespace Minuit2
 
 } // end namespace ROOT

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

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