Skip to content

Commit

Permalink
[RF] Avoid unneeded cached projection integrals in RooAddPdf
Browse files Browse the repository at this point in the history
By already normalizing the pdf components integrals to the right fit
range, one can elide some potentially expensive numeric correction
factor integrals.

The code that calculates the predicted number of events was updated
accordingly, since it containts also some integrals that depend on the
components normalization range.

Closes #12645.

A unit test that validates that there are no superfluous integrals in
the case of the GitHub reproducer was also added.
  • Loading branch information
guitargeek committed Jan 26, 2025
1 parent a5df5f1 commit c981070
Show file tree
Hide file tree
Showing 5 changed files with 97 additions and 9 deletions.
1 change: 1 addition & 0 deletions roofit/roofitcore/inc/RooFit/Detail/RooNormalizedPdf.h
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ class RooNormalizedPdf : public RooAbsPdf {
auto name = std::string(pdf.GetName()) + "_over_" + _normIntegral->GetName();
SetName(name.c_str());
SetTitle(name.c_str());
_normRange = pdf.normRange(); // so that e.g. RooAddPdf can query over what we are normalized
}

RooNormalizedPdf(const RooNormalizedPdf &other, const char *name)
Expand Down
11 changes: 7 additions & 4 deletions roofit/roofitcore/src/RooAddHelpers.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -56,14 +56,17 @@ AddCacheElem::AddCacheElem(RooAbsPdf const &addPdf, RooArgList const &pdfList, R
nset2.add(refCoefNormSet);
}

bool hasPdfWithCustomRange = false;
bool hasPdfWithDifferentRange = false;

// Fill with dummy unit RRVs for now
for (std::size_t i = 0; i < pdfList.size(); ++i) {
auto pdf = static_cast<const RooAbsPdf *>(pdfList.at(i));
auto coef = static_cast<const RooAbsReal *>(coefList.at(i));

hasPdfWithCustomRange |= pdf->normRange() != nullptr;
const std::string normRangeComponent = pdf->normRange() ? pdf->normRange() : "";
const bool componentHasDifferentNormRange = normRangeComponent != normRange;

hasPdfWithDifferentRange |= componentHasDifferentNormRange;

// Start with full list of dependents
RooArgSet supNSet(fullDepList);
Expand All @@ -88,7 +91,7 @@ AddCacheElem::AddCacheElem(RooAbsPdf const &addPdf, RooArgList const &pdfList, R
<< pdf->GetName() << std::endl;
}

if (!normRange.empty()) {
if (componentHasDifferentNormRange) {
auto snormTerm = std::unique_ptr<RooAbsReal>(pdf->createIntegral(nset2, nset2, normRange.c_str()));
if (snorm) {
auto oldSnorm = std::move(snorm);
Expand All @@ -109,7 +112,7 @@ AddCacheElem::AddCacheElem(RooAbsPdf const &addPdf, RooArgList const &pdfList, R

// *** PART 2 : Create projection coefficients ***

const bool projectCoefsForRangeReasons = !refCoefNormRange.empty() || !normRange.empty() || hasPdfWithCustomRange;
const bool projectCoefsForRangeReasons = !refCoefNormRange.empty() || !normRange.empty() || hasPdfWithDifferentRange;

// If no projections required stop here
if (refCoefNormSet.empty() && !projectCoefsForRangeReasons) {
Expand Down
50 changes: 46 additions & 4 deletions roofit/roofitcore/src/RooAddPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ An (enforced) condition for this assumption is that each \f$ \mathrm{PDF}_i \f$
#include <RooDataSet.h>
#include <RooGenericPdf.h>
#include <RooGlobalFunc.h>
#include <RooProdPdf.h>
#include <RooProduct.h>
#include <RooRatio.h>
#include <RooRealConstant.h>
Expand Down Expand Up @@ -823,7 +824,10 @@ std::unique_ptr<RooAbsReal> RooAddPdf::createExpectedEventsFunc(const RooArgSet
// too and we can't add them all in one go as owned objects of the
// final integral sum.
for (auto *pdf : static_range_cast<RooAbsPdf *>(_pdfList)) {
auto next = std::unique_ptr<RooAbsReal>{pdf->createIntegral(*nset, *nset, _normRange)};
auto integrl = std::unique_ptr<RooAbsReal>{pdf->createIntegral(*nset, *nset)};
auto formulaName = std::string(pdf->GetName()) + "_formulaVar";
auto next = std::make_unique<RooFormulaVar>(formulaName.c_str(), "1./x[0]", RooArgList{*integrl});
next->addOwnedComponents(std::move(integrl));
terms.add(*next);
if (owner)
next->addOwnedComponents(std::move(owner));
Expand Down Expand Up @@ -966,15 +970,53 @@ RooAddPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileCo
auto newArg = std::unique_ptr<RooAbsReal>{static_cast<RooAbsReal *>(Clone())};
ctx.markAsCompiled(*newArg);

// If we set the normalization ranges of the component pdfs to the
// normalization range of the RooAddPdf, the RooAddPdf doesn't need to
// compute as many correction factor integrals internally.
// Remember the previous normalizations ranges to reset later.
class ResetNormRangesRAII {
public:
ResetNormRangesRAII(RooAbsCollection const &pdfs, std::string const &normRange)
{
_componentPdfs.reserve(pdfs.size());
_oldNormRanges.reserve(pdfs.size());

bool isMultiRange = normRange.find(',') != std::string::npos;

for (auto *componentPdf : static_range_cast<RooAbsPdf *>(pdfs)) {

// The RooProdPdf is not able to deal with a multi-range _normRange
// for now, so skip the changing the range in this case.
bool changeRange = !(isMultiRange && dynamic_cast<RooProdPdf *>(componentPdf));

const char *old = componentPdf->normRange();
const char *newVal = changeRange ? normRange.c_str() : old;
componentPdf->setNormRange(newVal);
_componentPdfs.emplace_back(componentPdf);
_oldNormRanges.emplace_back(old ? old : "");
}
}
~ResetNormRangesRAII()
{
for (std::size_t i = 0; i < _componentPdfs.size(); ++i) {
_componentPdfs[i]->setNormRange(_oldNormRanges[i].c_str());
}
}

private:
std::vector<RooAbsPdf *> _componentPdfs;
std::vector<std::string> _oldNormRanges;
} resetNormRangesRAII{_pdfList, _normRange.Data()};

// In case conditional observables, e.g. p(x|y), the _refCoefNorm is set to
// all observables (x, y) and the normSet doesn't contain the conditional
// observables (so it only contains x in this example).

// If the _refCoefNorm is empty or it's equal to normSet anyway, this is not
// a conditional pdf and we don't need to do any transformation.
if(_refCoefNorm.empty() || normSet.equals(_refCoefNorm)) {
ctx.compileServers(*newArg, normSet);
return newArg;
if (_refCoefNorm.empty() || normSet.equals(_refCoefNorm)) {
ctx.compileServers(*newArg, normSet);
return newArg;
}

// In the conditional case, things become more complicated. The original
Expand Down
2 changes: 1 addition & 1 deletion roofit/roofitcore/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
# @author Patrick Bos, NL eScience Center, 2018

ROOT_ADD_GTEST(testSimple testSimple.cxx LIBRARIES RooFitCore)
ROOT_ADD_GTEST(testRooAddPdf testRooAddPdf.cxx LIBRARIES RooFitCore RooStats)
ROOT_ADD_GTEST(testRooAddPdf testRooAddPdf.cxx LIBRARIES RooFitCore RooFit RooStats)
ROOT_ADD_GTEST(testRooCacheManager testRooCacheManager.cxx LIBRARIES RooFitCore)
ROOT_ADD_GTEST(testRooCategory testRooCategory.cxx LIBRARIES RooFitCore)
ROOT_ADD_GTEST(testWorkspace testWorkspace.cxx LIBRARIES RooFitCore RooStats)
Expand Down
42 changes: 42 additions & 0 deletions roofit/roofitcore/test/testRooAddPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,14 @@
#include <RooConstVar.h>
#include <RooDataHist.h>
#include <RooFit/Evaluator.h>
#include <RooGaussian.h>
#include <RooGenericPdf.h>
#include <RooHelpers.h>
#include <RooHistPdf.h>
#include <RooMsgService.h>
#include <RooProdPdf.h>
#include <RooRealIntegral.h>
#include <RooUniform.h>
#include <RooWorkspace.h>

#include <RooStats/SPlot.h>
Expand Down Expand Up @@ -382,3 +385,42 @@ TEST(RooAddPdf, ImplicitDimensions)
EXPECT_DOUBLE_EQ(pdf.getVal(normSet), refVal);
EXPECT_DOUBLE_EQ(evaluator.run()[0], refVal);
}

// Make sure that there are no superfluous integrals when one does a ranged
// with equivalent coefficient reference range.
// Covers GitHub issue 12645.
TEST(RooAddPdf, IntegralsForRangedFitWithIdenticalCoefRange)
{
RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);

using namespace RooFit;

RooRealVar x("x", "", 0, 1);

RooGaussian g("g", "", x, 0.5, 0.2);
RooUniform u("u", "", x);

RooRealVar f("f", "", 0.5, 0, 1);
RooAddPdf a("a", "", {g, u}, {f});

std::unique_ptr<RooAbsData> dt{a.generate(x, 1000)};

x.setRange("limited", 0.2, 0.8);

std::unique_ptr<RooAbsReal> nll{a.createNLL(*dt, Range("limited"), SumCoefRange("limited"), EvalBackend("cpu"))};

RooArgList nodes;
nll->treeNodeServerList(&nodes);

int iIntegrals = 0;

for (auto *arg : nodes) {
if (dynamic_cast<RooRealIntegral const *>(arg)) {
++iIntegrals;
}
}

// We expect only two integral objects: one normalization integral for each
// of the component pdfs.
EXPECT_EQ(iIntegrals, 2);
}

0 comments on commit c981070

Please sign in to comment.