Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 2 additions & 0 deletions roofit/codegen/inc/RooFit/CodegenImpl.h
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ class RooRealSumFunc;
class RooRealSumPdf;
class RooRealVar;
class RooRecursiveFraction;
class RooSimultaneous;
class RooUniform;
class RooWrapperPdf;

Expand Down Expand Up @@ -126,6 +127,7 @@ void codegenImpl(RooRealSumFunc &arg, CodegenContext &ctx);
void codegenImpl(RooRealSumPdf &arg, CodegenContext &ctx);
void codegenImpl(RooRealVar &arg, CodegenContext &ctx);
void codegenImpl(RooRecursiveFraction &arg, CodegenContext &ctx);
void codegenImpl(RooSimultaneous &arg, CodegenContext &ctx);
void codegenImpl(RooStats::HistFactory::FlexibleInterpVar &arg, CodegenContext &ctx);
void codegenImpl(RooUniform &arg, CodegenContext &ctx);
void codegenImpl(RooWrapperPdf &arg, CodegenContext &ctx);
Expand Down
44 changes: 44 additions & 0 deletions roofit/codegen/src/CodegenImpl.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@
#include <RooRealSumPdf.h>
#include <RooRealVar.h>
#include <RooRecursiveFraction.h>
#include <RooSimultaneous.h>
#include <RooStats/HistFactory/FlexibleInterpVar.h>
#include <RooStats/HistFactory/ParamHistFunc.h>
#include <RooStats/HistFactory/PiecewiseInterpolation.h>
Expand Down Expand Up @@ -291,6 +292,49 @@ void codegenImpl(RooMultiPdf &arg, CodegenContext &ctx)
}
}

void codegenImpl(RooSimultaneous &arg, CodegenContext &ctx)
{
// A RooSimultaneous appears as a node in the compute graph only when its
// index category is not an observable (see the "switch mode" in
// RooSimultaneous::compileForNormSet()): its value is then the one of the
// component selected by the current category state, like for RooMultiPdf.
// With an observable index category, the likelihood is decomposed into
// per-channel terms and the RooSimultaneous itself is never translated.
if (arg.canBeExtended()) {
std::stringstream errorMsg;
errorMsg << "RooSimultaneous \"" << arg.GetName()
<< "\" with extendable components and a non-observable index category can't be translated, because the "
"scalar evaluation applies a relative yield weight that code generation does not implement yet.";
oocoutE(&arg, Minimization) << errorMsg.str() << std::endl;
throw std::runtime_error(errorMsg.str());
}
if (!arg.indexCat().isFundamental()) {
std::stringstream errorMsg;
errorMsg << "RooSimultaneous \"" << arg.GetName()
<< "\" with a derived, non-observable index category can't be translated.";
oocoutE(&arg, Minimization) << errorMsg.str() << std::endl;
throw std::runtime_error(errorMsg.str());
}

std::string const &indexExpr = ctx.getResult(arg.indexCat());

// Nested ternary expression that selects the pdf matching the category
// state, keyed on the state index numbers.
std::string expr;
std::size_t numStates = 0;
for (auto const &nameIdx : arg.indexCat()) {
RooAbsPdf *pdf = arg.getPdf(nameIdx.first.c_str());
if (!pdf) {
continue;
}
expr += "(" + indexExpr + " == " + std::to_string(nameIdx.second) + " ? (" + ctx.getResult(*pdf) + ") : ";
++numStates;
}
expr += "0.0" + std::string(numStates, ')');

ctx.addResult(&arg, expr);
}

// RooCategory index added.
void codegenImpl(RooCategory &arg, CodegenContext &ctx)
{
Expand Down
2 changes: 2 additions & 0 deletions roofit/roofitcore/inc/RooSimultaneous.h
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,8 @@ class RooSimultaneous : public RooAbsPdf {

std::unique_ptr<RooAbsArg> compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override;

bool indexCatIsObservable(RooArgSet const &vars) const;

protected:

void selectNormalization(const RooArgSet* depSet=nullptr, bool force=false) override ;
Expand Down
8 changes: 8 additions & 0 deletions roofit/roofitcore/src/BatchModeDataHelpers.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -236,6 +236,14 @@ RooFit::BatchModeDataHelpers::getDataSpans(RooAbsData const &data, std::string c
RooSimultaneous const *simPdf, bool skipZeroWeights,
bool takeGlobalObservablesFromData, std::stack<std::vector<double>> &buffers)
{
// If the index category is not among the data columns, the RooSimultaneous
// doesn't split the data into channels but acts as a "switch" selecting the
// component given by the current index state (analogous to RooMultiPdf).
// The data is then loaded like for any ordinary pdf.
if (simPdf && !simPdf->indexCatIsObservable(*data.get())) {
simPdf = nullptr;
}

std::vector<std::pair<std::string, RooAbsData const *>> datasets;
std::vector<bool> isBinnedL;
bool splitRange = false;
Expand Down
15 changes: 12 additions & 3 deletions roofit/roofitcore/src/FitHelpers.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -489,7 +489,12 @@ std::unique_ptr<RooAbsReal> createNLLNew(RooAbsPdf &pdf, RooAbsData &data, std::
RooAbsPdf &finalPdf = applyIntegrateBinsWrapping(pdf, data, integrateOverBinsPrecision, binSamplingPdfs);

RooArgList nllTerms;
if (auto *simPdf = dynamic_cast<RooSimultaneous *>(&finalPdf)) {
auto *simPdf = dynamic_cast<RooSimultaneous *>(&finalPdf);
// A RooSimultaneous whose index category is not among the data columns is
// a "switch" pdf selecting the component given by the current index state
// (analogous to RooMultiPdf): there are no channels to split the NLL into,
// so it is treated like an ordinary pdf.
if (simPdf && simPdf->indexCatIsObservable(*data.get())) {
nllTerms.addOwned(createSimultaneousNLL(*simPdf, isExtended, rangeName, offset));
} else {
RooNLLVarNew::Config cfg;
Expand Down Expand Up @@ -831,7 +836,8 @@ std::unique_ptr<RooAbsReal> createNLL(RooAbsPdf &pdf, RooAbsData &data, const Ro
RooArgSet normSet;
pdf.getObservables(data.get(), normSet);

if (dynamic_cast<RooSimultaneous const *>(&pdf)) {
auto *simPdfForProjDeps = dynamic_cast<RooSimultaneous const *>(&pdf);
if (simPdfForProjDeps && simPdfForProjDeps->indexCatIsObservable(normSet)) {
for (auto i : projDeps) {
auto res = normSet.find(i->GetName());
if (res != nullptr) {
Expand Down Expand Up @@ -1088,7 +1094,10 @@ std::unique_ptr<RooAbsReal> createChi2(RooAbsReal &real, RooDataHist &data, cons
applyIntegrateBinsWrapping(*pdfClone, data, pc.getDouble("integrate_bins"), binSamplingPdfs);

std::unique_ptr<RooAbsReal> chi2;
if (auto *simPdfClone = dynamic_cast<RooSimultaneous *>(&finalPdf)) {
auto *simPdfClone = dynamic_cast<RooSimultaneous *>(&finalPdf);
// Like in createNLLNew(): a "switch"-mode RooSimultaneous (index
// category not among the data columns) is treated as an ordinary pdf.
if (simPdfClone && simPdfClone->indexCatIsObservable(*data.get())) {
chi2 = std::unique_ptr<RooAbsReal>{dynamic_cast<RooAbsReal *>(
createSimultaneousChi2(*simPdfClone, rangeName ? rangeName : "", extended, etype).release())};
} else {
Expand Down
24 changes: 20 additions & 4 deletions roofit/roofitcore/src/RooAbsTestStatistic.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,21 @@ using std::endl, std::ostream;
/// ```
/// if the categories are called "pi0" and "gamma".

namespace {

/// A RooSimultaneous implies a simultaneous fit over the states of its index
/// category only if the index category is among the data columns. Otherwise,
/// it acts as a "switch" pdf that evaluates to the component selected by the
/// current index state (analogous to RooMultiPdf), and the test statistic has
/// to treat it like any ordinary pdf instead of splitting the data.
bool isSimultaneousFit(RooAbsReal &real, RooAbsData &data)
{
auto *simPdf = dynamic_cast<RooSimultaneous *>(&real);
return simPdf && simPdf->indexCatIsObservable(*data.get());
}

} // namespace

RooAbsTestStatistic::RooAbsTestStatistic(const char *name, const char *title, RooAbsReal& real, RooAbsData& data,
const RooArgSet& projDeps, RooAbsTestStatistic::Configuration const& cfg) :
RooAbsReal(name,title),
Expand All @@ -97,8 +112,8 @@ RooAbsTestStatistic::RooAbsTestStatistic(const char *name, const char *title, Ro
_addCoefRangeName(cfg.addCoefRangeName),
_splitRange(cfg.splitCutRange),
_verbose(cfg.verbose),
// Determine if RooAbsReal is a RooSimultaneous
_gofOpMode{(cfg.nCPU>1 || cfg.nCPU==-1) ? MPMaster : (dynamic_cast<RooSimultaneous*>(_func) ? SimMaster : Slave)},
// Determine if RooAbsReal implies a simultaneous fit over channels
_gofOpMode{(cfg.nCPU>1 || cfg.nCPU==-1) ? MPMaster : (isSimultaneousFit(real, data) ? SimMaster : Slave)},
_nEvents{data.numEntries()},
_nCPU(cfg.nCPU != -1 ? cfg.nCPU : 1),
_mpinterl(cfg.interleave),
Expand All @@ -123,8 +138,9 @@ RooAbsTestStatistic::RooAbsTestStatistic(const RooAbsTestStatistic& other, const
_addCoefRangeName(other._addCoefRangeName),
_splitRange(other._splitRange),
_verbose(other._verbose),
// Determine if RooAbsReal is a RooSimultaneous
_gofOpMode{(other._nCPU>1 || other._nCPU==-1) ? MPMaster : (dynamic_cast<RooSimultaneous*>(_func) ? SimMaster : Slave)},
// Determine if RooAbsReal implies a simultaneous fit over channels
_gofOpMode{(other._nCPU>1 || other._nCPU==-1) ? MPMaster
: (isSimultaneousFit(*other._func, *other._data) ? SimMaster : Slave)},
_nEvents{_data->numEntries()},
_nCPU(other._nCPU != -1 ? other._nCPU : 1),
_mpinterl(other._mpinterl),
Expand Down
27 changes: 27 additions & 0 deletions roofit/roofitcore/src/RooSimultaneous.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1250,6 +1250,22 @@ RooArgSet const& RooSimultaneous::flattenedCatList() const
return *_indexCatSet;
}

////////////////////////////////////////////////////////////////////////////////
/// Check if the index category is among the variables `vars`, matching by
/// name. For a RooSuperCategory index, any of its input categories counts.
///
/// If the index category is not among the observables of a fit, this
/// RooSimultaneous does not split the data into channels: it acts as a
/// "switch" that evaluates to the component selected by the current index
/// state, analogous to RooMultiPdf. Fitting infrastructure uses this check to
/// decide between the two modes.
bool RooSimultaneous::indexCatIsObservable(RooArgSet const &vars) const
{
RooArgSet catsAmongVars;
vars.selectCommon(flattenedCatList(), catsAmongVars);
return !catsAmongVars.empty();
}

namespace {

void markObs(RooAbsArg *arg, std::string const &prefix, RooArgSet const &normSet)
Expand Down Expand Up @@ -1284,6 +1300,17 @@ void prefixArgs(RooAbsArg *arg, std::string const &prefix, RooArgSet const &norm
std::unique_ptr<RooAbsArg>
RooSimultaneous::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const
{
if (!indexCatIsObservable(normSet)) {
// The index category is not an observable here: the RooSimultaneous
// acts as a plain "switch" that evaluates to the component selected by
// the current index state, analogous to RooMultiPdf. The channel
// observables are then the same as the ones of this pdf, so the
// channel-splitting compilation below (which renames the per-channel
// observables so they can be filled from split datasets) must not be
// used. Compile like an ordinary self-normalized pdf instead.
return RooAbsPdf::compileForNormSet(normSet, ctx);
}

std::unique_ptr<RooSimultaneous> newSimPdf{static_cast<RooSimultaneous *>(this->Clone())};

const char *rangeName = this->getStringAttribute("RangeName");
Expand Down
105 changes: 105 additions & 0 deletions roofit/roofitcore/test/testRooSimultaneous.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -829,6 +829,63 @@ TEST(RooSimultaneous, AsymmetryPlot)
}
}

/// A RooSimultaneous whose index category is not among the observables acts
/// as a "switch" that evaluates to the component selected by the current
/// index state, analogous to RooMultiPdf (see the discussion in GitHub issue
/// #22916). The batch backends used to silently compute wrong values for such
/// a pdf nested inside another model, because the compilation for the
/// likelihood assumed that the dataset can be split by the index category.
TEST(RooSimultaneous, ParameterIndexSwitchMode)
{
using namespace RooFit;

RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
RooRandom::randomGenerator()->SetSeed(1337ul);

RooRealVar x("x", "x", 0, 50);
RooRealVar lam("lam", "lam", -0.04, -0.1, -0.01);
RooExponential expo("expo", "expo", x, lam);
RooRealVar c0("c0", "c0", -0.5, -1.0, 1.0);
RooRealVar c1("c1", "c1", 0.2, -1.0, 1.0);
RooChebychev cheb("cheb", "cheb", x, {c0, c1});

RooRealVar mean("mean", "mean", 25, 20, 30);
RooRealVar sigma("sigma", "sigma", 2, 0.5, 5);
RooGaussian gauss("gauss", "gauss", x, mean, sigma);
RooRealVar frac("frac", "frac", 0.8, 0.0, 1.0);

RooCategory cat("cat", "cat", {{"expo", 0}, {"cheb", 1}});
RooSimultaneous sim("sim", "sim", cat);
sim.addPdf(expo, "expo");
sim.addPdf(cheb, "cheb");
RooAddPdf model("model", "model", {sim, gauss}, frac);

// Reference models with the switch resolved by hand.
RooAddPdf refModel0("refModel0", "refModel0", {expo, gauss}, frac);
RooAddPdf refModel1("refModel1", "refModel1", {cheb, gauss}, frac);

std::unique_ptr<RooDataSet> data{refModel0.generate(x, 500)};

std::vector<RooFit::EvalBackend> backends;
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
backends.push_back(RooFit::EvalBackend::Legacy());
#endif
backends.push_back(RooFit::EvalBackend::Cpu());
backends.push_back(RooFit::EvalBackend::CodegenNoGrad());

for (auto &backend : backends) {
cat.setIndex(0);
std::unique_ptr<RooAbsReal> nll{model.createNLL(*data, backend)};
std::unique_ptr<RooAbsReal> refNll0{refModel0.createNLL(*data, backend)};
std::unique_ptr<RooAbsReal> refNll1{refModel1.createNLL(*data, backend)};

cat.setIndex(0);
EXPECT_THAT(nll->getVal(), RelativeNear(refNll0->getVal(), 1e-10)) << backend.name() << ", index 0";
cat.setIndex(1);
EXPECT_THAT(nll->getVal(), RelativeNear(refNll1->getVal(), 1e-10)) << backend.name() << ", index 1";
}
}

/// Regression test for a value-server link corruption: the index category is
/// a value server of the RooSimultaneous via the index category proxy, and
/// fixAddCoefNormalization() additionally stores it in the non-propagating
Expand Down Expand Up @@ -858,3 +915,51 @@ TEST(RooSimultaneous, RepeatedFixAddCoefNormalization)
EXPECT_TRUE(observables.find(sample)) << "after fixAddCoefNormalization() call " << i + 1;
}
}

/// A "switch"-mode RooSimultaneous (see the ParameterIndexSwitchMode test
/// above) must also work as the top-level pdf of a fit: creating the NLL used
/// to fail on all backends, because the fitting infrastructure unconditionally
/// tried to split the dataset by the index category. Now the index category is
/// treated as a parameter of the fit if it is not among the data columns,
/// making the RooSimultaneous usable in place of RooMultiPdf for the discrete
/// profiling method.
TEST(RooSimultaneous, ParameterIndexTopLevelNLL)
{
using namespace RooFit;

RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
RooRandom::randomGenerator()->SetSeed(1337ul);

RooRealVar x("x", "x", 0, 50);
RooRealVar lam("lam", "lam", -0.04, -0.1, -0.01);
RooExponential expo("expo", "expo", x, lam);
RooRealVar c0("c0", "c0", -0.5, -1.0, 1.0);
RooRealVar c1("c1", "c1", 0.2, -1.0, 1.0);
RooChebychev cheb("cheb", "cheb", x, {c0, c1});

RooCategory cat("cat", "cat", {{"expo", 0}, {"cheb", 1}});
RooSimultaneous sim("sim", "sim", cat);
sim.addPdf(expo, "expo");
sim.addPdf(cheb, "cheb");

std::unique_ptr<RooDataSet> data{expo.generate(x, 500)};

std::vector<RooFit::EvalBackend> backends;
#ifdef ROOFIT_LEGACY_EVAL_BACKEND
backends.push_back(RooFit::EvalBackend::Legacy());
#endif
backends.push_back(RooFit::EvalBackend::Cpu());
backends.push_back(RooFit::EvalBackend::CodegenNoGrad());

for (auto &backend : backends) {
cat.setIndex(0);
std::unique_ptr<RooAbsReal> nll{sim.createNLL(*data, backend)};
std::unique_ptr<RooAbsReal> refNll0{expo.createNLL(*data, backend)};
std::unique_ptr<RooAbsReal> refNll1{cheb.createNLL(*data, backend)};

cat.setIndex(0);
EXPECT_THAT(nll->getVal(), RelativeNear(refNll0->getVal(), 1e-10)) << backend.name() << ", index 0";
cat.setIndex(1);
EXPECT_THAT(nll->getVal(), RelativeNear(refNll1->getVal(), 1e-10)) << backend.name() << ", index 1";
}
}
Loading