diff --git a/roofit/codegen/inc/RooFit/CodegenImpl.h b/roofit/codegen/inc/RooFit/CodegenImpl.h index 33bf102cef6f8..9e02b392cf17a 100644 --- a/roofit/codegen/inc/RooFit/CodegenImpl.h +++ b/roofit/codegen/inc/RooFit/CodegenImpl.h @@ -60,6 +60,7 @@ class RooRealSumFunc; class RooRealSumPdf; class RooRealVar; class RooRecursiveFraction; +class RooSimultaneous; class RooUniform; class RooWrapperPdf; @@ -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); diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index b1f3545352318..bc3bdc3556d3e 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -55,6 +55,7 @@ #include #include #include +#include #include #include #include @@ -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) { diff --git a/roofit/roofitcore/inc/RooSimultaneous.h b/roofit/roofitcore/inc/RooSimultaneous.h index f4c51531fa627..b29ac80a3df34 100644 --- a/roofit/roofitcore/inc/RooSimultaneous.h +++ b/roofit/roofitcore/inc/RooSimultaneous.h @@ -107,6 +107,8 @@ class RooSimultaneous : public RooAbsPdf { std::unique_ptr 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 ; diff --git a/roofit/roofitcore/src/BatchModeDataHelpers.cxx b/roofit/roofitcore/src/BatchModeDataHelpers.cxx index d40b03fcbb879..eb6aeeff74c7e 100644 --- a/roofit/roofitcore/src/BatchModeDataHelpers.cxx +++ b/roofit/roofitcore/src/BatchModeDataHelpers.cxx @@ -236,6 +236,14 @@ RooFit::BatchModeDataHelpers::getDataSpans(RooAbsData const &data, std::string c RooSimultaneous const *simPdf, bool skipZeroWeights, bool takeGlobalObservablesFromData, std::stack> &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> datasets; std::vector isBinnedL; bool splitRange = false; diff --git a/roofit/roofitcore/src/FitHelpers.cxx b/roofit/roofitcore/src/FitHelpers.cxx index 04bbc9724e75c..f1a33a1159d47 100644 --- a/roofit/roofitcore/src/FitHelpers.cxx +++ b/roofit/roofitcore/src/FitHelpers.cxx @@ -489,7 +489,12 @@ std::unique_ptr createNLLNew(RooAbsPdf &pdf, RooAbsData &data, std:: RooAbsPdf &finalPdf = applyIntegrateBinsWrapping(pdf, data, integrateOverBinsPrecision, binSamplingPdfs); RooArgList nllTerms; - if (auto *simPdf = dynamic_cast(&finalPdf)) { + auto *simPdf = dynamic_cast(&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; @@ -831,7 +836,8 @@ std::unique_ptr createNLL(RooAbsPdf &pdf, RooAbsData &data, const Ro RooArgSet normSet; pdf.getObservables(data.get(), normSet); - if (dynamic_cast(&pdf)) { + auto *simPdfForProjDeps = dynamic_cast(&pdf); + if (simPdfForProjDeps && simPdfForProjDeps->indexCatIsObservable(normSet)) { for (auto i : projDeps) { auto res = normSet.find(i->GetName()); if (res != nullptr) { @@ -1088,7 +1094,10 @@ std::unique_ptr createChi2(RooAbsReal &real, RooDataHist &data, cons applyIntegrateBinsWrapping(*pdfClone, data, pc.getDouble("integrate_bins"), binSamplingPdfs); std::unique_ptr chi2; - if (auto *simPdfClone = dynamic_cast(&finalPdf)) { + auto *simPdfClone = dynamic_cast(&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{dynamic_cast( createSimultaneousChi2(*simPdfClone, rangeName ? rangeName : "", extended, etype).release())}; } else { diff --git a/roofit/roofitcore/src/RooAbsTestStatistic.cxx b/roofit/roofitcore/src/RooAbsTestStatistic.cxx index b2dbe20e6f89d..ff5a6549b09c7 100644 --- a/roofit/roofitcore/src/RooAbsTestStatistic.cxx +++ b/roofit/roofitcore/src/RooAbsTestStatistic.cxx @@ -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(&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), @@ -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(_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), @@ -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(_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), diff --git a/roofit/roofitcore/src/RooSimultaneous.cxx b/roofit/roofitcore/src/RooSimultaneous.cxx index 56b7503dee643..9e4257000075a 100644 --- a/roofit/roofitcore/src/RooSimultaneous.cxx +++ b/roofit/roofitcore/src/RooSimultaneous.cxx @@ -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) @@ -1284,6 +1300,17 @@ void prefixArgs(RooAbsArg *arg, std::string const &prefix, RooArgSet const &norm std::unique_ptr 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 newSimPdf{static_cast(this->Clone())}; const char *rangeName = this->getStringAttribute("RangeName"); diff --git a/roofit/roofitcore/test/testRooSimultaneous.cxx b/roofit/roofitcore/test/testRooSimultaneous.cxx index dee568994e40f..6792bd866fdf2 100644 --- a/roofit/roofitcore/test/testRooSimultaneous.cxx +++ b/roofit/roofitcore/test/testRooSimultaneous.cxx @@ -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 data{refModel0.generate(x, 500)}; + + std::vector 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 nll{model.createNLL(*data, backend)}; + std::unique_ptr refNll0{refModel0.createNLL(*data, backend)}; + std::unique_ptr 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 @@ -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 data{expo.generate(x, 500)}; + + std::vector 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 nll{sim.createNLL(*data, backend)}; + std::unique_ptr refNll0{expo.createNLL(*data, backend)}; + std::unique_ptr 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"; + } +}