[Minuit2] Merge FCNAdapter with -GradAdapter and use only std::function

The base classes `FCNBase` and `FCNGradientBase` were already unified in
71806b3, and the Minuit code can be further simplified by merging
also the derived adapter class for `std::function` objects.

This will also make it easier for the users, because then don't have to
use different classes depending on whether they pass external gradients
or not.

We also make the new `FCNAdapter` purely `std::function` based, instead
of using a template parameter. This makes is possible to use Minuit2
directly from Python, as we can benefit from the `std::function`
pythonization.

It also makes the Minuit2 code more consistent, because `std::function`
was already used for the G2 (second derivatives) and the Hessian.

Author: guitargeek

math/minuit2/CMakeLists.txt

@@ -26,7 +26,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT)
       Minuit2/ExternalInternalGradientCalculator.h
       Minuit2/FCNAdapter.h
       Minuit2/FCNBase.h
-      Minuit2/FCNGradAdapter.h
       Minuit2/FCNGradientBase.h
       Minuit2/FumiliBuilder.h
       Minuit2/FumiliChi2FCN.h

math/minuit2/inc/Minuit2/FCNAdapter.h

@@ -15,6 +15,7 @@
 #include <ROOT/RSpan.hxx>
 
 #include <vector>
+#include <functional>
 
 namespace ROOT {
 
@@ -23,29 +24,94 @@ namespace Minuit2 {
 /**
 
 
-template wrapped class for adapting to FCNBase signature
+Wrapper class for std::function objects, adapting to FCNBase signature.
 
 @author Lorenzo Moneta
 
 @ingroup Minuit
 
 */
 
-template <class Function>
 class FCNAdapter : public FCNBase {
 
 public:
-   FCNAdapter(const Function &f, double up = 1.) : fFunc(f), fUp(up) {}
+   FCNAdapter(std::function<double(double const *)> f, double up = 1.) : fUp(up), fFunc(std::move(f)) {}
+
+   bool HasGradient() const override { return bool(fGradFunc); }
+   bool HasG2() const override { return bool(fG2Func); }
+   bool HasHessian() const override { return bool(fHessianFunc); }
+
+   double operator()(std::vector<double> const &v) const override { return fFunc(v.data()); }
 
-   double operator()(std::vector<double> const& v) const override { return fFunc.operator()(&v[0]); }
-   double operator()(const double *v) const { return fFunc.operator()(v); }
    double Up() const override { return fUp; }
 
+   std::vector<double> Gradient(std::vector<double> const &v) const override
+   {
+      std::vector<double> output(v.size());
+      fGradFunc(v.data(), output.data());
+      return output;
+   }
+
+   /// return second derivatives (diagonal of the Hessian matrix)
+   std::vector<double> G2(std::vector<double> const &x) const override
+   {
+      std::vector<double> output;
+      if (fG2Func)
+         return fG2Func(x);
+      if (fHessianFunc) {
+         std::size_t n = x.size();
+         output.resize(n);
+         if (fHessian.empty())
+            fHessian.resize(n * n);
+         fHessianFunc(x, fHessian.data());
+         if (!fHessian.empty()) {
+            // get diagonal element of h
+            for (unsigned int i = 0; i < n; i++)
+               output[i] = fHessian[i * n + i];
+         }
+      }
+      return output;
+   }
+
+   /// compute Hessian. Return Hessian as a std::vector of size(n*n)
+   std::vector<double> Hessian(std::vector<double> const &x) const override
+   {
+      std::vector<double> output;
+      if (fHessianFunc) {
+         std::size_t n = x.size();
+         output.resize(n * n);
+         bool ret = fHessianFunc(x, output.data());
+         if (!ret) {
+            output.clear();
+            fHessianFunc = nullptr;
+         }
+      }
+
+      return output;
+   }
+
+   void SetGradientFunction(std::function<void(double const *, double *)> f) { fGradFunc = std::move(f); }
+   void SetG2Function(std::function<std::vector<double>(std::vector<double> const &)> f) { fG2Func = std::move(f); }
+   void SetHessianFunction(std::function<bool(std::vector<double> const &, double *)> f)
+   {
+      fHessianFunc = std::move(f);
+   }
+
    void SetErrorDef(double up) override { fUp = up; }
 
 private:
-   const Function &fFunc;
-   double fUp;
+   using Function = std::function<double(double const *)>;
+   using GradFunction = std::function<void(double const *, double *)>;
+   using G2Function = std::function<std::vector<double>(std::vector<double> const &)>;
+   using HessianFunction = std::function<bool(std::vector<double> const &, double *)>;
+
+   double fUp = 1.;
+   mutable std::vector<double> fHessian;
+
+   Function fFunc;
+   GradFunction fGradFunc;
+   G2Function fG2Func;
+   mutable HessianFunction fHessianFunc;
 };
 
 } // end namespace Minuit2

math/minuit2/inc/Minuit2/FCNBase.h

@@ -113,6 +113,10 @@ class FCNBase {
    virtual bool HasGradient() const { return false; }
 
    virtual std::vector<double> Gradient(std::vector<double> const&) const { return {}; }
+
+   /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
    virtual std::vector<double> GradientWithPrevResult(std::vector<double> const& parameters, double * /*previous_grad*/,
                                                       double * /*previous_g2*/, double * /*previous_gstep*/) const
    {

math/minuit2/inc/Minuit2/FCNGradAdapter.h

@@ -16,7 +16,6 @@ set(MINUIT2_HEADERS
     ExternalInternalGradientCalculator.h
     FCNAdapter.h
     FCNBase.h
-    FCNGradAdapter.h
     FCNGradientBase.h
     FumiliBuilder.h
     FumiliChi2FCN.h

math/minuit2/src/CMakeLists.txt

@@ -19,7 +19,6 @@
 
 #include "Minuit2/FCNAdapter.h"
 #include "Minuit2/FumiliFCNAdapter.h"
-#include "Minuit2/FCNGradAdapter.h"
 #include "Minuit2/FunctionMinimum.h"
 #include "Minuit2/MnMigrad.h"
 #include "Minuit2/MnMinos.h"
@@ -375,11 +374,14 @@ void Minuit2Minimizer::SetFunction(const ROOT::Math::IMultiGenFunction &func)
    fDim = func.NDim();
    const bool hasGrad = func.HasGradient();
    if (!fUseFumili) {
-      if (hasGrad)
-         fMinuitFCN = std::make_unique<ROOT::Minuit2::FCNGradAdapter<ROOT::Math::IMultiGradFunction>>(
-            dynamic_cast<ROOT::Math::IMultiGradFunction const &>(func), ErrorDef());
-      else
-         fMinuitFCN = std::make_unique<ROOT::Minuit2::FCNAdapter<ROOT::Math::IMultiGenFunction>>(func, ErrorDef());
+      auto lambdaFunc = [&func](double const *params) { return func(params); };
+      auto adapter = std::make_unique<ROOT::Minuit2::FCNAdapter>(lambdaFunc, ErrorDef());
+      if (hasGrad) {
+         auto const &gradFunc = dynamic_cast<ROOT::Math::IMultiGradFunction const &>(func);
+         auto lambdaGrad = [&gradFunc](double const *params, double *grad) { return gradFunc.Gradient(params, grad); };
+         adapter->SetGradientFunction(lambdaGrad);
+      }
+      fMinuitFCN = std::move(adapter);
       return;
    }
    if (hasGrad) {
@@ -409,7 +411,7 @@ void Minuit2Minimizer::SetHessianFunction(std::function<bool(std::span<const dou
 {
    // for Fumili not supported for the time being
    if (fUseFumili) return;
-   auto fcn = static_cast<ROOT::Minuit2::FCNGradAdapter<ROOT::Math::IMultiGradFunction> *>(fMinuitFCN.get());
+   auto fcn = static_cast<ROOT::Minuit2::FCNAdapter *>(fMinuitFCN.get());
    if (!fcn) return;
    fcn->SetHessianFunction(hfunc);
 }