[Minuit2] Performance and thread safety improvements plus refactoring

math/minuit2/CMakeLists.txt

@@ -40,7 +40,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT)
       Minuit2/FumiliStandardMaximumLikelihoodFCN.h
       Minuit2/FunctionGradient.h
       Minuit2/FunctionMinimum.h
-      Minuit2/GenericFunction.h
       Minuit2/GradientCalculator.h
       Minuit2/HessianGradientCalculator.h
       Minuit2/InitialGradientCalculator.h
@@ -150,7 +149,6 @@ if(CMAKE_PROJECT_NAME STREQUAL ROOT)
       src/MnStrategy.cxx
       src/MnTiny.cxx
       src/MnTraceObject.cxx
-      src/MnUserFcn.cxx
       src/MnUserParameterState.cxx
       src/MnUserParameters.cxx
       src/MnUserTransformation.cxx

math/minuit2/inc/LinkDef.h

@@ -13,9 +13,6 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
-//#pragma link C++ namespace ROOT::Minuit2;
-
-#pragma link C++ class ROOT::Minuit2::GenericFunction;
 #pragma link C++ class ROOT::Minuit2::FCNBase;
 #pragma link C++ class ROOT::Minuit2::FCNGradientBase;
 #pragma link C++ class ROOT::Minuit2::FumiliFCNBase;

math/minuit2/inc/Minuit2/FCNBase.h

@@ -12,8 +12,6 @@
 
 #include "Minuit2/MnConfig.h"
 
-#include "Minuit2/GenericFunction.h"
-
 #include <ROOT/RSpan.hxx>
 
 #include <vector>
@@ -48,10 +46,12 @@ Interface (abstract class) defining the function to be minimized, which has to b
 
  */
 
-class FCNBase : public GenericFunction {
+class FCNBase {
 
 public:
 
+   virtual ~FCNBase() = default;
+
    /**
 
       The meaning of the vector of parameters is of course defined by the user,
@@ -74,7 +74,7 @@ class FCNBase : public GenericFunction {
 
    */
 
-   double operator()(std::vector<double> const &v) const override = 0;
+   virtual double operator()(std::vector<double> const &v) const = 0;
 
    /**
 

math/minuit2/inc/Minuit2/FumiliBuilder.h

@@ -13,7 +13,6 @@
 #include "Minuit2/MinimumBuilder.h"
 #include "Minuit2/VariableMetricEDMEstimator.h"
 #include "Minuit2/FumiliErrorUpdator.h"
-#include "Minuit2/MnFcn.h"
 #include "Minuit2/FunctionMinimum.h"
 
 #include <vector>
@@ -22,6 +21,8 @@ namespace ROOT {
 
 namespace Minuit2 {
 
+class MnFcn;
+
 /**
 
 Builds the FunctionMinimum using the Fumili method.

math/minuit2/inc/Minuit2/MnFcn.h

@@ -10,13 +10,17 @@
 #ifndef ROOT_Minuit2_MnFcn
 #define ROOT_Minuit2_MnFcn
 
-#include "Minuit2/MnConfig.h"
+#include "Minuit2/FCNBase.h"
 #include "Minuit2/MnMatrix.h"
 
+#include <vector>
+
 namespace ROOT {
 
 namespace Minuit2 {
 
+class MnUserTransformation;
+
 class FCNBase;
 /**
    Wrapper class to FCNBase interface used internally by Minuit.
@@ -30,27 +34,51 @@ class FCNBase;
 class MnFcn {
 
 public:
-   /// constructor of
    explicit MnFcn(const FCNBase &fcn, int ncall = 0) : fFCN(fcn), fNumCall(ncall) {}
+   explicit MnFcn(const FCNBase &fcn, const MnUserTransformation &trafo, int ncall = 0)
+      : fFCN(fcn), fNumCall(ncall), fTransform(&trafo)
+   {
+   }
 
-   virtual ~MnFcn();
-
-   virtual double operator()(const MnAlgebraicVector &) const;
    unsigned int NumOfCalls() const { return fNumCall; }
 
-   //
-   // forward interface
-   //
-   double ErrorDef() const;
-   double Up() const;
+   double ErrorDef() const
+   {
+      return fFCN.Up();
+   }
+
+   double Up() const
+   {
+      return fFCN.Up();
+   }
 
    const FCNBase &Fcn() const { return fFCN; }
 
+   // Access the parameter transformations.
+   // For internal use in the Minuit2 implementation.
+   const MnUserTransformation *Trafo() const { return fTransform; }
+
+   double CallWithTransformedParams(std::vector<double> const &vpar) const;
+   double CallWithoutDoingTrafo(const MnAlgebraicVector &) const;
+
 private:
    const FCNBase &fFCN;
-
-protected:
    mutable int fNumCall;
+   const MnUserTransformation *fTransform = nullptr;
+};
+
+// Helper class to call the MnFcn, caching the transformed parameters if necessary.
+class MnFcnCaller {
+public:
+   MnFcnCaller(const MnFcn &mfcn);
+
+   double operator()(const MnAlgebraicVector &v);
+
+private:
+   MnFcn const &fMfcn;
+   bool fDoInt2ext = false;
+   std::vector<double> fLastInput;
+   std::vector<double> fVpar;
 };
 
 } // namespace Minuit2

math/minuit2/inc/Minuit2/MnUserFcn.h

@@ -1,51 +1,14 @@
-// @(#)root/minuit2:$Id$
-// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
-
-/**********************************************************************
- *                                                                    *
- * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
- *                                                                    *
- **********************************************************************/
-
 #ifndef ROOT_Minuit2_MnUserFcn
 #define ROOT_Minuit2_MnUserFcn
 
-#include "Minuit2/MnFcn.h"
-
-#include <vector>
+#include <Minuit2/MnFcn.h>
 
 namespace ROOT {
-
 namespace Minuit2 {
 
-class MnUserTransformation;
-
-/**
-   Wrapper used by Minuit of FCN interface
-   containing a reference to the transformation object
- */
-class MnUserFcn : public MnFcn {
-
-public:
-   MnUserFcn(const FCNBase &fcn, const MnUserTransformation &trafo, int ncall = 0)
-      : MnFcn(fcn, ncall), fTransform(trafo)
-   {
-   }
-
-   double operator()(const MnAlgebraicVector &) const override;
-
-   // Access the parameter transformations.
-   // For internal use in the Minuit2 implementation.
-   const MnUserTransformation &transform() const { return fTransform; }
-
-   double callWithTransformedParams(std::vector<double> const &vpar) const;
-
-private:
-   const MnUserTransformation &fTransform;
-};
-
-} // namespace Minuit2
+using MnUserFcn = MnFcn;
 
+}
 } // namespace ROOT
 
-#endif // ROOT_Minuit2_MnUserFcn
+#endif

math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h

@@ -14,17 +14,12 @@
 
 #include "Minuit2/GradientCalculator.h"
 
-#include <ROOT/RSpan.hxx>
-
-#include <vector>
-
 namespace ROOT {
 
 namespace Minuit2 {
 
 class MnFcn;
 class MnUserTransformation;
-class MnMachinePrecision;
 class MnStrategy;
 
 /**
@@ -43,15 +38,6 @@ class Numerical2PGradientCalculator : public GradientCalculator {
 
    FunctionGradient operator()(const MinimumParameters &, const FunctionGradient &) const override;
 
-   const MnFcn &Fcn() const { return fFcn; }
-   const MnUserTransformation &Trafo() const { return fTransformation; }
-   const MnMachinePrecision &Precision() const;
-   const MnStrategy &Strategy() const { return fStrategy; }
-
-   unsigned int Ncycle() const;
-   double StepTolerance() const;
-   double GradTolerance() const;
-
 private:
    const MnFcn &fFcn;
    const MnUserTransformation &fTransformation;

math/minuit2/src/CMakeLists.txt

@@ -30,7 +30,6 @@ set(MINUIT2_HEADERS
     FumiliStandardMaximumLikelihoodFCN.h
     FunctionGradient.h
     FunctionMinimum.h
-    GenericFunction.h
     GradientCalculator.h
     HessianGradientCalculator.h
     InitialGradientCalculator.h
@@ -76,7 +75,6 @@ set(MINUIT2_HEADERS
     MnTiny.h
     MnTraceObject.h
     MnUserCovariance.h
-    MnUserFcn.h
     MnUserParameterState.h
     MnUserParameters.h
     MnUserTransformation.h
@@ -139,7 +137,6 @@ set(MINUIT2_SOURCES
     MnStrategy.cxx
     MnTiny.cxx
     MnTraceObject.cxx
-    MnUserFcn.cxx
     MnUserParameterState.cxx
     MnUserParameters.cxx
     MnUserTransformation.cxx

math/minuit2/src/CombinedMinimumBuilder.cxx

@@ -7,7 +7,6 @@
  *                                                                    *
  **********************************************************************/
 
-#include "Minuit2/AnalyticalGradientCalculator.h"
 #include "Minuit2/CombinedMinimumBuilder.h"
 #include "Minuit2/FunctionMinimum.h"
 #include "Minuit2/MnStrategy.h"
@@ -37,10 +36,7 @@ FunctionMinimum CombinedMinimumBuilder::Minimum(const MnFcn &fcn, const Gradient
 
          return min1;
       }
-      // check if gradient calculator is analytical
-      auto agc = dynamic_cast<const AnalyticalGradientCalculator *>(&gc);
-      MinimumSeed seed1 = (agc) ?
-         fVMMinimizer.SeedGenerator()(fcn, *agc, min1.UserState(), str) : fVMMinimizer.SeedGenerator()(fcn, gc, min1.UserState(), str);
+      MinimumSeed seed1 = fVMMinimizer.SeedGenerator()(fcn, gc, min1.UserState(), str);
 
       FunctionMinimum min2 = fVMMinimizer.Builder().Minimum(fcn, gc, seed1, str, maxfcn, edmval);
       if (!min2.IsValid()) {

math/minuit2/src/FumiliBuilder.cxx

@@ -10,7 +10,6 @@
 #include "Minuit2/FumiliBuilder.h"
 #include "Minuit2/FumiliStandardMaximumLikelihoodFCN.h"
 #include "Minuit2/GradientCalculator.h"
-//#include "Minuit2/Numerical2PGradientCalculator.h"
 #include "Minuit2/MinimumState.h"
 #include "Minuit2/MinimumError.h"
 #include "Minuit2/FunctionGradient.h"
@@ -237,6 +236,8 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
 
    double lambda = (doLineSearch) ? 0.001 : 0;
 
+   MnFcnCaller fcnCaller{fcn};
+
    do {
 
       //     const MinimumState& s0 = result.back();
@@ -269,7 +270,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
       // take a full step
 
       //evaluate function only if doing a line search
-      double fval2 = (doLineSearch) ? fcn(s0.Vec() + step) : 0;
+      double fval2 = (doLineSearch) ? fcnCaller(s0.Vec() + step) : 0;
       MinimumParameters p(s0.Vec() + step, fval2);
 
       // check that taking the full step does not deteriorate minimum
@@ -324,7 +325,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
 
          // if the proposed point (newton step) is inside the trust region radius accept it
          if (norm <= delta) {
-            p = MinimumParameters(s0.Vec() + step, fcn(s0.Vec() + step));
+            p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
             print.Debug("Accept full Newton step - it is inside TR ",delta);
          } else {
             //step = - (delta/norm) * step;
@@ -403,7 +404,7 @@ FunctionMinimum FumiliBuilder::Minimum(const MnFcn &fcn, const GradientCalculato
             }
             print.Debug("New accepted step is ",step);
 
-            p = MinimumParameters(s0.Vec() + step, fcn(s0.Vec() + step));
+            p = MinimumParameters(s0.Vec() + step, fcnCaller(s0.Vec() + step));
             norm = delta;
             gdel = inner_product(step, s0.Gradient().Grad());
          }