[RF] Test RooFit Hessians with Clad

Author: guitargeek

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