Skip to content

[RF] Replace legacy backend reference in tests with analytical results #22092

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
guitargeek wants to merge 1 commit into
base: master
Choose a base branch
from
+223 −127
Open

[RF] Replace legacy backend reference in tests with analytical results #22092

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ ValAndError getValAndError(RooArgSet const &parsFinal, const char *name)
return {var.getVal(), var.getError()};
};

#ifndef ROOFIT_LEGACY_EVAL_BACKEND
std::vector<double> getParamVals(RooAbsMinimizerFcn &fcn)
{
std::vector<double> values(fcn.getNDim());
Expand All @@ -73,6 +74,7 @@ std::vector<double> getParamVals(RooAbsMinimizerFcn &fcn)

return values;
}
#endif

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

#ifdef ROOFIT_LEGACY_EVAL_BACKEND_
TEST_P(LikelihoodGradientJobTest, GaussianND)
{
// do a minimization, but now using GradMinimizer and its MP version
Expand Down Expand Up @@ -315,6 +318,7 @@ TEST_P(LikelihoodGradientJobTest, GaussianND)
EXPECT_EQ(std0[ix], std1[ix]);
}
}
#endif

INSTANTIATE_TEST_SUITE_P(NworkersSeed, LikelihoodGradientJobTest,
::testing::Combine(::testing::Values(1, 2, 3), // number of workers
Expand Down Expand Up @@ -579,6 +583,7 @@ TEST_P(LikelihoodGradientJobTest, Gaussian1DAlsoWithLikelihoodJob)
}
#undef EXPECT_NEAR_REL

#ifndef ROOFIT_LEGACY_EVAL_BACKEND
class LikelihoodGradientJobErrorTest
: public ::testing::TestWithParam<std::tuple<std::size_t, std::size_t, bool, bool>> {
void SetUp() override
Expand Down Expand Up @@ -967,3 +972,4 @@ TEST(MinuitFcnGrad, DISABLED_CompareToRooMinimizerFcn)
RFMP::Config::LikelihoodJob::defaultNEventTasks = RFMP::Config::LikelihoodJob::automaticNEventTasks;
RFMP::Config::LikelihoodJob::defaultNComponentTasks = RFMP::Config::LikelihoodJob::automaticNComponentTasks;
}
#endif
55 changes: 39 additions & 16 deletions roofit/roofitcore/test/testRooAbsPdf.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,9 +96,6 @@ TEST_P(FitTest, AsymptoticallyCorrectErrors)
// evaluated in batch mode and data size is greater than one, the batch mode
// will inform that a batched evaluation function is missing.
//
// This test is disabled if the legacy backend is not available, because then
// we don't have any reference to compare to.
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
TEST(RooAbsPdf, ConditionalFitBatchMode)
{
using namespace RooFit;
Expand Down Expand Up @@ -132,41 +129,67 @@ TEST(RooAbsPdf, ConditionalFitBatchMode)

auto data = makeFakeDataXY();

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

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

std::vector<bool> expectFastEvaluationsWarnings{true, false};
double sumX = 0.0;
double sumY = 0.0;
for (int i = 0; i < data->numEntries(); ++i) {
const RooArgSet *row = data->get(i);
sumX += row->getRealValue("x");
sumY += row->getRealValue("y");
}
const double nEntries = data->numEntries();

// For each event, the conditional log-likelihood term is
// log Poisson(x_i; factor * y_i) = x_i log(factor*y_i) - factor*y_i + const
// (the Poisson normalisation integral over x is unity by construction of
// the wide x range). Setting d(NLL)/d(factor) = 0 gives:
// model1: factor_MLE = sum_i x_i / sum_i y_i (mean depends on y)
// model2: factor_MLE = sum_i x_i / N (mean is just factor)
// and the standard error from the inverse Hessian is in both cases
// sigma(factor) = sqrt(sum_i x_i) / (sum_i y_i or N).
const std::vector<double> expectedFactor{sumX / sumY, sumX / nEntries};
const std::vector<double> expectedFactorErr{std::sqrt(sumX) / sumY, std::sqrt(sumX) / nEntries};
const std::vector<bool> expectFastEvaluationsWarnings{true, false};

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

std::vector<std::unique_ptr<RooFitResult>> fitResults;

RooHelpers::HijackMessageStream hijack(RooFit::INFO, RooFit::FastEvaluations);

for (auto evalBackend : {EvalBackend::Legacy(), EvalBackend::Cpu()}) {
factor.setVal(1.0);
factor.setError(0.0);
fitResults.emplace_back(model->fitTo(*data, ConditionalObservables(y), Save(), PrintLevel(-1), evalBackend));
if (verbose) {
fitResults.back()->Print();
}
factor.setVal(1.0);
factor.setError(0.0);
auto fitResult = std::unique_ptr<RooFitResult>{
model->fitTo(*data, ConditionalObservables(y), Save(), PrintLevel(-1), EvalBackend::Cpu())};
if (verbose) {
fitResult->Print();
}

EXPECT_TRUE(fitResults[1]->isIdentical(*fitResults[0]));
auto *factorFinal = static_cast<RooRealVar *>(fitResult->floatParsFinal().find("factor"));
ASSERT_NE(factorFinal, nullptr);
EXPECT_NEAR(factorFinal->getVal(), expectedFactor[iMean], 1e-4 * expectedFactor[iMean])
<< "value mismatch for " << model->GetName();
EXPECT_NEAR(factorFinal->getError(), expectedFactorErr[iMean], 1e-3 * expectedFactorErr[iMean])
<< "error mismatch for " << model->GetName();

EXPECT_EQ(hijack.str().find("does not implement the faster batch") != std::string::npos,
expectFastEvaluationsWarnings[iMean])
<< "Stream contents: " << hijack.str();
++iMean;
}
}
#endif

// ROOT-9530: RooFit side-band fit inconsistent with fit to full range
TEST_P(FitTest, MultiRangeFit)
Expand Down
47 changes: 33 additions & 14 deletions roofit/roofitcore/test/testRooSimultaneous.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@

#include "gtest_wrapper.h"

#include <cmath>
#include <memory>

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

const double nBkgA_nom = 9000;
const double nBkgB_nom = 10000;

RooRealVar x("x", "Observable", 100, 150);
x.setRange("fitRange", 100, 130);
const double cA = -0.06;
const double cB = -0.09;

const double xMin = 100;
const double xMax = 150;
const double xFitMax = 130;

RooRealVar x("x", "Observable", xMin, xMax);
x.setRange("fitRange", xMin, xFitMax);

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

RooExponential expA("expA", "", x, RooFit::RooConst(-0.06));
RooExponential expA("expA", "", x, RooFit::RooConst(cA));
RooAddPdf modelA("modelA", "", {expA}, {nBkgA});

RooExponential expB("expB", "", x, RooFit::RooConst(-0.09));
RooExponential expB("expB", "", x, RooFit::RooConst(cB));
RooAddPdf modelB("modelB", "", {expB}, {nBkgB});

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

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

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

RooArgSet params;
RooArgSet paramsSnap;
simPdf.getParameters(combData->get(), params);
params.snapshot(paramsSnap);
ASSERT_NE(fitResult, nullptr);
EXPECT_EQ(fitResult->status(), 0);

auto integ = [](double c, double a, double b) { return (std::exp(c * b) - std::exp(c * a)) / c; };

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

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

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

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

/// GitHub issue #20383.
Expand Down
Loading
Loading