lmoneta · lmoneta · Feb 28, 2022
diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt
@@ -139,6 +139,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
     RooMappedCategory.h
     RooMCIntegrator.h
     RooMCStudy.h
+    RooCladMinimizerFcn.h
     RooMinimizerFcn.h
     RooMinimizer.h
     RooMoment.h
@@ -368,6 +369,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
     src/RooMCIntegrator.cxx
     src/RooMCStudy.cxx
     src/RooMinimizer.cxx
+    src/RooCladMinimizerFcn.cxx
     src/RooMinimizerFcn.cxx
     src/RooMoment.cxx
     src/RooMPSentinel.cxx
@@ -460,6 +462,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
     src/TestStatistics/LikelihoodGradientWrapper.cxx
     src/TestStatistics/LikelihoodWrapper.cxx
     src/TestStatistics/LikelihoodSerial.cxx
+    src/TestStatistics/LikelihoodGradientClad.cxx
     src/TestStatistics/MinuitFcnGrad.cxx
     src/TestStatistics/RooAbsL.cxx
     src/TestStatistics/RooBinnedL.cxx

@@ -374,6 +374,8 @@ class RooAbsReal : public RooAbsArg {
   static void setHideOffset(Bool_t flag);
   static Bool_t hideOffset() ;
 
+  virtual void evaluateGradient(double*) const {}
+
 protected:
   // Hook for objects with normalization-dependent parameters interpretation
   virtual void selectNormalization(const RooArgSet* depSet=0, Bool_t force=kFALSE) ;

@@ -0,0 +1,71 @@
+/*****************************************************************************
+ * Project: RooFit                                                           *
+ * Package: RooFitCore                                                       *
+ * @(#)root/roofitcore:$Id$
+ * Authors:                                                                  *
+ *   AL, Alfio Lazzaro,   INFN Milan,        [email protected]        *
+ *   PB, Patrick Bos, Netherlands eScience Center, [email protected]   *
+ *                                                                           *
+ *                                                                           *
+ * Redistribution and use in source and binary forms,                        *
+ * with or without modification, are permitted according to the terms        *
+ * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
+ *****************************************************************************/
+
+#ifndef roofit_roofitcore_RooCladMinimizerFcn_h
+#define roofit_roofitcore_RooCladMinimizerFcn_h
+
+#include "Math/IFunction.h"
+#include "Fit/ParameterSettings.h"
+#include "Fit/FitResult.h"
+
+#include "RooAbsReal.h"
+#include "RooArgList.h"
+
+#include <fstream>
+#include <vector>
+
+#include <RooAbsMinimizerFcn.h>
+
+template<typename T> class TMatrixTSym;
+using TMatrixDSym = TMatrixTSym<double>;
+
+// forward declaration
+class RooMinimizer;
+
+// class RooCladMinimizerFcn : public RooAbsMinimizerFcn, public ROOT::Math::IBaseFunctionMultiDim {
+class RooCladMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbsMinimizerFcn {
+
+public:
+   RooCladMinimizerFcn(RooAbsReal *funct, RooMinimizer *context, bool verbose = false);
+   RooCladMinimizerFcn(const RooCladMinimizerFcn &other);
+   ~RooCladMinimizerFcn() override;
+
+   ROOT::Math::IBaseFunctionMultiDim *Clone() const override;
+
+   /// IMultiGradFunction overrides necessary for Minuit
+   void Gradient(const double *x, double *grad) const override;
+   void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2,
+                               double *previous_gstep) const override;
+
+   unsigned int NDim() const override { return getNDim(); }
+
+   std::string getFunctionName() const override;
+   std::string getFunctionTitle() const override;
+
+   void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override;
+
+   void setOffsetting(Bool_t flag) override;
+
+private:
+   /// IMultiGradFunction overrides necessary for Minuit
+   double DoDerivative(const double *x, unsigned int icoord) const override;
+   double DoDerivativeWithPrevResult(const double *x, unsigned int i_component, double *previous_grad,
+                                     double *previous_g2, double *previous_gstep) const override;
+
+   double DoEval(const double *x) const override;
+
+   RooAbsReal *_funct;
+};
+
+#endif
@@ -30,7 +30,7 @@ namespace TestStatistics {
 class RooAbsL;
 struct WrapperCalculationCleanFlags;
 
-enum class LikelihoodGradientMode { multiprocess };
+enum class LikelihoodGradientMode { multiprocess, clad };
 
 class LikelihoodGradientWrapper {
 public:

@@ -95,6 +95,8 @@ class RooAbsL {
    virtual ROOT::Math::KahanSum<double>
    evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) = 0;
 
+   virtual void evaluateGradient(double*) const {}
+
    // necessary from MinuitFcnGrad to reach likelihood properties:
    virtual RooArgSet *getParameters();
 

@@ -48,7 +48,7 @@ class RooPlot ;
 
 class RooMinimizer : public TObject {
 public:
-  enum class FcnMode { classic, gradient, generic_wrapper };
+  enum class FcnMode { classic, gradient, generic_wrapper, clad };
 
   explicit RooMinimizer(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic);
   explicit RooMinimizer(std::shared_ptr<RooFit::TestStatistics::RooAbsL> likelihood,

@@ -0,0 +1,165 @@
+/*****************************************************************************
+ * Project: RooFit                                                           *
+ * Package: RooFitCore                                                       *
+ * @(#)root/roofitcore:$Id$
+ * Authors:                                                                  *
+ *   AL, Alfio Lazzaro,   INFN Milan,        [email protected]        *
+ *   PB, Patrick Bos, Netherlands eScience Center, [email protected]   *
+ *                                                                           *
+ *                                                                           *
+ * Redistribution and use in source and binary forms,                        *
+ * with or without modification, are permitted according to the terms        *
+ * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
+ *****************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+/// \class RooCladMinimizerFcn
+/// RooCladMinimizerFcn is an interface to the ROOT::Math::IBaseFunctionMultiDim,
+/// a function that ROOT's minimisers use to carry out minimisations.
+///
+
+#include "RooCladMinimizerFcn.h"
+
+#include "RooAbsArg.h"
+#include "RooAbsPdf.h"
+#include "RooArgSet.h"
+#include "RooRealVar.h"
+#include "RooAbsRealLValue.h"
+#include "RooMsgService.h"
+#include "RooMinimizer.h"
+#include "RooNaNPacker.h"
+
+#include "TClass.h"
+#include "TMatrixDSym.h"
+
+#include <fstream>
+#include <iomanip>
+
+using namespace std;
+
+
+namespace {
+
+// Helper function that wraps RooAbsArg::getParameters and directly returns the
+// output RooArgSet. To be used in the initializer list of the RooCladMinimizerFcn
+// constructor.
+RooArgSet getParameters(RooAbsReal const& funct) {
+    RooArgSet out;
+    funct.getParameters(nullptr, out);
+    return out;
+}
+
+} // namespace
+
+
+RooCladMinimizerFcn::RooCladMinimizerFcn(RooAbsReal *funct, RooMinimizer* context,
+            bool verbose) :
+  RooAbsMinimizerFcn(getParameters(*funct), context, verbose), _funct(funct)
+{}
+
+
+
+RooCladMinimizerFcn::RooCladMinimizerFcn(const RooCladMinimizerFcn& other) : ROOT::Math::IMultiGradFunction(other), RooAbsMinimizerFcn(other),
+  _funct(other._funct)
+{}
+
+
+RooCladMinimizerFcn::~RooCladMinimizerFcn()
+{}
+
+
+ROOT::Math::IBaseFunctionMultiDim* RooCladMinimizerFcn::Clone() const
+{
+  return new RooCladMinimizerFcn(*this) ;
+}
+
+void RooCladMinimizerFcn::setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt)
+{
+   _funct->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt);
+}
+
+/// Evaluate function given the parameters in `x`.
+double RooCladMinimizerFcn::DoEval(const double *x) const {
+
+  // Set the parameter values for this iteration
+  for (unsigned index = 0; index < _nDim; index++) {
+    if (_logfile) (*_logfile) << x[index] << " " ;
+    SetPdfParamVal(index,x[index]);
+  }
+
+  // Calculate the function for these parameters
+  RooAbsReal::setHideOffset(kFALSE) ;
+  double fvalue = _funct->getVal();
+  RooAbsReal::setHideOffset(kTRUE) ;
+
+  if (!std::isfinite(fvalue) || RooAbsReal::numEvalErrors() > 0 || fvalue > 1e30) {
+    printEvalErrors();
+    RooAbsReal::clearEvalErrorLog() ;
+    _numBadNLL++ ;
+
+    if (_doEvalErrorWall) {
+      const double badness = RooNaNPacker::unpackNaN(fvalue);
+      fvalue = (std::isfinite(_maxFCN) ? _maxFCN : 0.) + _recoverFromNaNStrength * badness;
+    }
+  } else {
+    if (_evalCounter > 0 && _evalCounter == _numBadNLL) {
+      // This is the first time we get a valid function value; while before, the
+      // function was always invalid. For invalid  cases, we returned values > 0.
+      // Now, we offset valid values such that they are < 0.
+      _funcOffset = -fvalue;
+    }
+    fvalue += _funcOffset;
+    _maxFCN = std::max(fvalue, _maxFCN);
+  }
+
+  // Optional logging
+  if (_logfile)
+    (*_logfile) << setprecision(15) << fvalue << setprecision(4) << endl;
+  if (_verbose) {
+    cout << "\nprevFCN" << (_funct->isOffsetting()?"-offset":"") << " = " << setprecision(10)
+         << fvalue << setprecision(4) << "  " ;
+    cout.flush() ;
+  }
+
+  _evalCounter++ ;
+
+  return fvalue;
+}
+
+std::string RooCladMinimizerFcn::getFunctionName() const
+{
+   return _funct->GetName();
+}
+
+std::string RooCladMinimizerFcn::getFunctionTitle() const
+{
+   return _funct->GetTitle();
+}
+
+void RooCladMinimizerFcn::setOffsetting(Bool_t flag)
+{
+   _funct->enableOffsetting(flag);
+}
+
+void RooCladMinimizerFcn::Gradient(const double *x, double *grad) const {
+   DoEval(x);
+   _funct->evaluateGradient(grad);
+}
+
+void RooCladMinimizerFcn::GradientWithPrevResult(const double *x, double *grad, double * /*previous_grad*/, double * /*previous_g2*/,
+                            double * /*previous_gstep*/) const {
+   Gradient(x, grad);
+}
+
+double RooCladMinimizerFcn::DoDerivative(const double *x, unsigned int icoord) const {
+   double grad[_nDim];
+   Gradient(x, grad);
+   return grad[icoord];
+}
+
+double RooCladMinimizerFcn::DoDerivativeWithPrevResult(const double *x, unsigned int i_component, double *previous_grad,
+                                  double *previous_g2, double *previous_gstep) const {
+   double grad[_nDim];
+   Gradient(x, grad);
+   return grad[i_component];
+}
@@ -51,6 +51,7 @@ automatic PDF optimization.
 #include "RooMsgService.h"
 #include "RooPlot.h"
 #include "RooMinimizerFcn.h"
+#include "RooCladMinimizerFcn.h"
 #include "RooGradMinimizerFcn.h"
 #include "RooFitResult.h"
 #include "TestStatistics/MinuitFcnGrad.h"
@@ -108,6 +109,10 @@ RooMinimizer::RooMinimizer(RooAbsReal &function, FcnMode fcnMode) : _fcnMode(fcn
       _fcn = new RooMinimizerFcn(&function, this, _verbose);
       break;
    }
+   case FcnMode::clad: {
+      _fcn = new RooCladMinimizerFcn(&function, this, _verbose);
+      break;
+   }
    case FcnMode::gradient: {
       _fcn = new RooGradMinimizerFcn(&function, this, _verbose);
       setMinimizerType("Minuit2");
@@ -297,6 +302,10 @@ bool RooMinimizer::fitFcn() const {
       ret = _theFitter->FitFCN(*dynamic_cast<RooMinimizerFcn *>(_fcn));
       break;
    }
+   case FcnMode::clad: {
+      ret = _theFitter->FitFCN(*dynamic_cast<RooCladMinimizerFcn *>(_fcn));
+      break;
+   }
    case FcnMode::gradient: {
       ret = _theFitter->FitFCN(*dynamic_cast<RooGradMinimizerFcn *>(_fcn));
       break;
@@ -835,6 +844,9 @@ ROOT::Math::IMultiGenFunction* RooMinimizer::getMultiGenFcn() const
       case FcnMode::classic: {
          return static_cast<ROOT::Math::IMultiGenFunction *>(dynamic_cast<RooMinimizerFcn *>(_fcn));
       }
+      case FcnMode::clad: {
+         return static_cast<ROOT::Math::IMultiGenFunction *>(dynamic_cast<RooCladMinimizerFcn *>(_fcn));
+      }
       case FcnMode::gradient: {
          return static_cast<ROOT::Math::IMultiGenFunction *>(dynamic_cast<RooGradMinimizerFcn *>(_fcn));
       }
@@ -859,6 +871,9 @@ const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const
       case FcnMode::gradient: {
          return static_cast<RooAbsMinimizerFcn *>(dynamic_cast<RooGradMinimizerFcn *>(getFitterMultiGenFcn()));
       }
+      case FcnMode::clad: {
+         return static_cast<RooAbsMinimizerFcn *>(dynamic_cast<RooCladMinimizerFcn *>(getFitterMultiGenFcn()));
+      }
       case FcnMode::generic_wrapper: {
          return static_cast<RooAbsMinimizerFcn *>(dynamic_cast<RooFit::TestStatistics::MinuitFcnGrad *>(getFitterMultiGenFcn()));
       }