[RF] Test RooFit Hessians with Clad

Author: guitargeek

interpreter/cling/tools/plugins/clad/CMakeLists.txt

@@ -74,7 +74,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.2)
+  list(APPEND _clad_extra_settings GIT_TAG master)
 endif()
 
 ## list(APPEND _clad_patches_list "patch1.patch" "patch2.patch")

roofit/histfactory/test/testHistFactory.cxx

@@ -9,17 +9,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>
@@ -618,9 +620,21 @@ TEST_P(HFFixtureFit, Fit)
       }
 
       using namespace RooFit;
-      std::unique_ptr<RooFitResult> fitResult{simPdf->fitTo(*data, evalBackend, Optimize(constTermOptimization),
-                                                            GlobalObservables(*mc->GetGlobalObservables()), Save(),
-                                                            PrintLevel(verbose ? 1 : -1))};
+      std::unique_ptr<RooAbsReal> nll{simPdf->createNLL(*data, evalBackend, Optimize(constTermOptimization),
+                                                        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");

roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h

@@ -180,7 +180,7 @@ double constraintSum(DoubleArray comp, unsigned int compSize)
 {
    double sum = 0;
 #if defined(__CLING__) && defined(R__HAS_CLAD)
-#pragma clad checkpoint loop
+//#pragma clad checkpoint loop
 #endif
    for (unsigned int i = 0; i < compSize; i++) {
       sum -= std::log(comp[i]);
@@ -370,7 +370,7 @@ double flexibleInterp(unsigned int code, ParamsArray params, unsigned int n, Dou
 {
    double total = nominal;
 #if defined(__CLING__) && defined(R__HAS_CLAD)
-#pragma clad checkpoint loop
+//#pragma clad checkpoint loop
 #endif
    for (std::size_t i = 0; i < n; ++i) {
       total += flexibleInterpSingle(code, low[i], high[i], boundary, nominal, params[i], total);

roofit/roofitcore/src/RooFit/CodegenContext.cxx

@@ -186,7 +186,7 @@ std::unique_ptr<CodegenContext::LoopScope> CodegenContext::beginLoop(RooAbsArg c
    }
 
    // Make sure that the name of this variable doesn't clash with other stuff
-   addToCodeBody(in, "#pragma clad checkpoint loop\n");
+   //addToCodeBody(in, "#pragma clad checkpoint loop\n");
    addToCodeBody(in, "for(int " + idx + " = 0; " + idx + " < obs[" + std::to_string(2 * firstObsIdx + 1) + "]; " + idx +
                         "++) {\n");
 

roofit/roofitcore/test/testRooFuncWrapper.cxx

@@ -129,14 +129,16 @@ class FactoryTest : public testing::TestWithParam<FactoryTestParams> {
    std::unique_ptr<RooHelpers::LocalChangeMsgLevel> _changeMsgLvl;
 };
 
-std::unique_ptr<RooFitResult> runMinimizer(RooAbsReal &absReal, bool useGradient = true)
+std::unique_ptr<RooFitResult> runMinimizer(RooAbsReal &absReal, bool useAD = true)
 {
    RooMinimizer::Config cfg;
-   cfg.useGradient = useGradient;
+   cfg.useGradient = useAD;
+   cfg.useHessian = useAD;
    RooMinimizer m{absReal, cfg};
    m.setPrintLevel(-1);
    m.setStrategy(0);
    m.minimize("Minuit2");
+   m.hesse(); // to use the analytical Hessian for error estimation
    return std::unique_ptr<RooFitResult>{m.save()};
 }
 
@@ -164,6 +166,9 @@ TEST_P(FactoryTest, NLLFit)
    // We want to use the generated code also for the nominal likelihood. Like
    // this, we make sure to validate also the NLL values of the generated code.
    static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nllFunc).setUseGeneratedFunctionCode(true);
+   // Let's also validate the Hessian
+   static_cast<RooFit::Experimental::RooEvaluatorWrapper&>(*nllFunc).generateHessian();
+   static_cast<RooFit::Experimental::RooEvaluatorWrapper&>(*nllFunc).writeDebugMacro(_params._name);
 
    double tol = _params._fitResultTolerance;
 
@@ -214,6 +219,9 @@ TEST_P(FactoryTest, NLLFit)
    // because for very small correlations it's usually not the same within the
    // relative tolerance because you would compare two small values that are
    // only different from zero because of noise.
+   // The error tolerance is also quite large, because the analytical Hessian
+   // gives numerically quite different results.
+   double errorTol = 0.1;
    EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, tol));
    EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, tol));
 }