Skip to content

Commit b2f7a03

Browse files
committed
[RF] Replace legacy backend reference in tests with analytical results
Some RooFit unit tests previously used the legacy evaluation backend to produce reference fit results, either as a direct comparison or guarded by `#ifdef ROOFIT_LEGACY_EVAL_BACKEND`. They are rewritten to compare against analytically-computed reference values, so they now exercise the CPU backend whether or not the legacy backend is built. Also, add missing `#ifdef ROOFIT_LEGACY_EVAL_BACKEND` guard to some tests in `testLikelihoodGradientJob.cxx`, because these tests check for bitwise equality with the legacy backend fit results.
1 parent f5c22fd commit b2f7a03

4 files changed

Lines changed: 223 additions & 127 deletions

File tree

roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@ ValAndError getValAndError(RooArgSet const &parsFinal, const char *name)
6363
return {var.getVal(), var.getError()};
6464
};
6565

66+
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
6667
std::vector<double> getParamVals(RooAbsMinimizerFcn &fcn)
6768
{
6869
std::vector<double> values(fcn.getNDim());
@@ -73,6 +74,7 @@ std::vector<double> getParamVals(RooAbsMinimizerFcn &fcn)
7374

7475
return values;
7576
}
77+
#endif
7678

7779
std::unique_ptr<RooFitResult> runMinimizer(RooAbsReal &nll, bool offsetting)
7880
{
@@ -235,6 +237,7 @@ TEST(LikelihoodGradientJob, RepeatMigrad)
235237
m1.minimize("Minuit2", "migrad");
236238
}
237239

240+
#ifdef ROOFIT_LEGACY_EVAL_BACKEND_
238241
TEST_P(LikelihoodGradientJobTest, GaussianND)
239242
{
240243
// do a minimization, but now using GradMinimizer and its MP version
@@ -315,6 +318,7 @@ TEST_P(LikelihoodGradientJobTest, GaussianND)
315318
EXPECT_EQ(std0[ix], std1[ix]);
316319
}
317320
}
321+
#endif
318322

319323
INSTANTIATE_TEST_SUITE_P(NworkersSeed, LikelihoodGradientJobTest,
320324
::testing::Combine(::testing::Values(1, 2, 3), // number of workers
@@ -579,6 +583,7 @@ TEST_P(LikelihoodGradientJobTest, Gaussian1DAlsoWithLikelihoodJob)
579583
}
580584
#undef EXPECT_NEAR_REL
581585

586+
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
582587
class LikelihoodGradientJobErrorTest
583588
: public ::testing::TestWithParam<std::tuple<std::size_t, std::size_t, bool, bool>> {
584589
void SetUp() override
@@ -967,3 +972,4 @@ TEST(MinuitFcnGrad, DISABLED_CompareToRooMinimizerFcn)
967972
RFMP::Config::LikelihoodJob::defaultNEventTasks = RFMP::Config::LikelihoodJob::automaticNEventTasks;
968973
RFMP::Config::LikelihoodJob::defaultNComponentTasks = RFMP::Config::LikelihoodJob::automaticNComponentTasks;
969974
}
975+
#endif

roofit/roofitcore/test/testRooAbsPdf.cxx

Lines changed: 39 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -96,9 +96,6 @@ TEST_P(FitTest, AsymptoticallyCorrectErrors)
9696
// evaluated in batch mode and data size is greater than one, the batch mode
9797
// will inform that a batched evaluation function is missing.
9898
//
99-
// This test is disabled if the legacy backend is not available, because then
100-
// we don't have any reference to compare to.
101-
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
10299
TEST(RooAbsPdf, ConditionalFitBatchMode)
103100
{
104101
using namespace RooFit;
@@ -132,41 +129,67 @@ TEST(RooAbsPdf, ConditionalFitBatchMode)
132129

133130
auto data = makeFakeDataXY();
134131

132+
// The model x range is wider than the support of the data so that the
133+
// Poisson normalisation integral over x is unity to high precision for any
134+
// mean value encountered. With that simplification, the conditional MLE for
135+
// `factor` has a closed form (see below) and the legacy evaluation backend
136+
// is no longer needed as a reference.
135137
RooWorkspace ws;
136138
ws.factory("Product::mean1({factor[1.0, 0.0, 10.0], y[1.0, 5]})");
137139
ws.factory("Product::mean2({factor})");
138-
ws.factory("Poisson::model1(x[0, 10], mean1)");
140+
ws.factory("Poisson::model1(x[0, 30], mean1)");
139141
ws.factory("Poisson::model2(x, mean2)");
140142

141143
RooRealVar &factor = *ws.var("factor");
142144
RooRealVar &y = *ws.var("y");
143145

144-
std::vector<bool> expectFastEvaluationsWarnings{true, false};
146+
double sumX = 0.0;
147+
double sumY = 0.0;
148+
for (int i = 0; i < data->numEntries(); ++i) {
149+
const RooArgSet *row = data->get(i);
150+
sumX += row->getRealValue("x");
151+
sumY += row->getRealValue("y");
152+
}
153+
const double nEntries = data->numEntries();
154+
155+
// For each event, the conditional log-likelihood term is
156+
// log Poisson(x_i; factor * y_i) = x_i log(factor*y_i) - factor*y_i + const
157+
// (the Poisson normalisation integral over x is unity by construction of
158+
// the wide x range). Setting d(NLL)/d(factor) = 0 gives:
159+
// model1: factor_MLE = sum_i x_i / sum_i y_i (mean depends on y)
160+
// model2: factor_MLE = sum_i x_i / N (mean is just factor)
161+
// and the standard error from the inverse Hessian is in both cases
162+
// sigma(factor) = sqrt(sum_i x_i) / (sum_i y_i or N).
163+
const std::vector<double> expectedFactor{sumX / sumY, sumX / nEntries};
164+
const std::vector<double> expectedFactorErr{std::sqrt(sumX) / sumY, std::sqrt(sumX) / nEntries};
165+
const std::vector<bool> expectFastEvaluationsWarnings{true, false};
145166

146167
int iMean = 0;
147168
for (RooAbsPdf *model : {ws.pdf("model1"), ws.pdf("model2")}) {
148169

149-
std::vector<std::unique_ptr<RooFitResult>> fitResults;
150-
151170
RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::FastEvaluations);
152171

153-
for (auto evalBackend : {EvalBackend::Legacy(), EvalBackend::Cpu()}) {
154-
factor.setVal(1.0);
155-
factor.setError(0.0);
156-
fitResults.emplace_back(model->fitTo(*data, ConditionalObservables(y), Save(), PrintLevel(-1), evalBackend));
157-
if (verbose) {
158-
fitResults.back()->Print();
159-
}
172+
factor.setVal(1.0);
173+
factor.setError(0.0);
174+
auto fitResult = std::unique_ptr<RooFitResult>{
175+
model->fitTo(*data, ConditionalObservables(y), Save(), PrintLevel(-1), EvalBackend::Cpu())};
176+
if (verbose) {
177+
fitResult->Print();
160178
}
161179

162-
EXPECT_TRUE(fitResults[1]->isIdentical(*fitResults[0]));
180+
auto *factorFinal = static_cast<RooRealVar *>(fitResult->floatParsFinal().find("factor"));
181+
ASSERT_NE(factorFinal, nullptr);
182+
EXPECT_NEAR(factorFinal->getVal(), expectedFactor[iMean], 1e-4 * expectedFactor[iMean])
183+
<< "value mismatch for " << model->GetName();
184+
EXPECT_NEAR(factorFinal->getError(), expectedFactorErr[iMean], 1e-3 * expectedFactorErr[iMean])
185+
<< "error mismatch for " << model->GetName();
186+
163187
EXPECT_EQ(hijack.str().find("does not implement the faster batch") != std::string::npos,
164188
expectFastEvaluationsWarnings[iMean])
165189
<< "Stream contents: " << hijack.str();
166190
++iMean;
167191
}
168192
}
169-
#endif
170193

171194
// ROOT-9530: RooFit side-band fit inconsistent with fit to full range
172195
TEST_P(FitTest, MultiRangeFit)

roofit/roofitcore/test/testRooSimultaneous.cxx

Lines changed: 33 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525

2626
#include "gtest_wrapper.h"
2727

28+
#include <cmath>
2829
#include <memory>
2930

3031
/// Forum issue
@@ -519,23 +520,36 @@ TEST_P(TestStatisticTest, RooSimultaneousSingleChannelCrossCheckWithCondVar)
519520
/// GitHub issue #18718.
520521
/// Make sure that we can do a ranged fit on an extended RooAddPdf in a
521522
/// RooSimultaneous with the new CPU backend.
523+
///
524+
/// The reference value is computed analytically: each channel has a single
525+
/// fixed-shape exponential, so the extended-MLE for the yield reduces to
526+
/// `N_obs_in_range * I_full / I_range`, where the integrals run over the full
527+
/// observable range and the fit range respectively. This avoids depending on
528+
/// the legacy evaluation backend, which is not always built.
522529
TEST(RooSimultaneous, RangedExtendedRooAddPdf)
523530
{
524531
RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
525532

526533
const double nBkgA_nom = 9000;
527534
const double nBkgB_nom = 10000;
528535

529-
RooRealVar x("x", "Observable", 100, 150);
530-
x.setRange("fitRange", 100, 130);
536+
const double cA = -0.06;
537+
const double cB = -0.09;
538+
539+
const double xMin = 100;
540+
const double xMax = 150;
541+
const double xFitMax = 130;
542+
543+
RooRealVar x("x", "Observable", xMin, xMax);
544+
x.setRange("fitRange", xMin, xFitMax);
531545

532546
RooRealVar nBkgA("nBkgA", "", nBkgA_nom, 0.8 * nBkgA_nom, 1.2 * nBkgA_nom);
533547
RooRealVar nBkgB("nBkgB", "", nBkgB_nom, 0.8 * nBkgB_nom, 1.2 * nBkgB_nom);
534548

535-
RooExponential expA("expA", "", x, RooFit::RooConst(-0.06));
549+
RooExponential expA("expA", "", x, RooFit::RooConst(cA));
536550
RooAddPdf modelA("modelA", "", {expA}, {nBkgA});
537551

538-
RooExponential expB("expB", "", x, RooFit::RooConst(-0.09));
552+
RooExponential expB("expB", "", x, RooFit::RooConst(cB));
539553
RooAddPdf modelB("modelB", "", {expB}, {nBkgB});
540554

541555
RooCategory runCat("runCat", "", {{"RunA", 0}, {"RunB", 1}});
@@ -546,21 +560,26 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
546560

547561
std::unique_ptr<RooDataSet> combData{simPdf.generate(RooArgSet(x, runCat), Extended())};
548562

549-
using Res = std::unique_ptr<RooFitResult>;
563+
std::unique_ptr<RooFitResult> fitResult{
564+
simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Cpu()), PrintLevel(-1))};
550565

551-
RooArgSet params;
552-
RooArgSet paramsSnap;
553-
simPdf.getParameters(combData->get(), params);
554-
params.snapshot(paramsSnap);
566+
ASSERT_NE(fitResult, nullptr);
567+
EXPECT_EQ(fitResult->status(), 0);
568+
569+
auto integ = [](double c, double a, double b) { return (std::exp(c * b) - std::exp(c * a)) / c; };
555570

556-
Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Cpu()), PrintLevel(-1))};
571+
const double ratioA = integ(cA, xMin, xMax) / integ(cA, xMin, xFitMax);
572+
const double ratioB = integ(cB, xMin, xMax) / integ(cB, xMin, xFitMax);
557573

558-
params.assign(paramsSnap);
574+
const double nBkgA_ref = combData->sumEntries("runCat==0", "fitRange") * ratioA;
575+
const double nBkgB_ref = combData->sumEntries("runCat==1", "fitRange") * ratioB;
559576

560-
Res fitResultRef{
561-
simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Legacy()), PrintLevel(-1))};
577+
auto getFinal = [&](const char *name) { return static_cast<RooRealVar *>(fitResult->floatParsFinal().find(name)); };
562578

563-
EXPECT_TRUE(fitResult->isIdentical(*fitResultRef));
579+
// Tolerance accounts for MINUIT's default convergence precision. The
580+
// fit-vs-analytical mismatch caused by the bug in #18718 was several percent.
581+
EXPECT_NEAR(getFinal("nBkgA")->getVal(), nBkgA_ref, 1e-3 * nBkgA_ref);
582+
EXPECT_NEAR(getFinal("nBkgB")->getVal(), nBkgB_ref, 1e-3 * nBkgB_ref);
564583
}
565584

566585
/// GitHub issue #20383.

0 commit comments

Comments
 (0)