From e2935d98efe66599056a861bd86392a3585665d5 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 20 Apr 2026 21:40:13 +0200 Subject: [PATCH 1/4] [RF] Use RooDataHist interfaces with index input arguments Use RooDataHist interfaces that take index input arguments, which is safer because they don't rely on the cached last index in the RooDataHist. Addresses one bullet point in https://github.com/root-project/root/issues/6557: "Correct interface of RooAbsData and derived classes to use e.g. std::size_t for indexing events. int doesn't make sense." As a followup, we could consider to deprecate the interfaces that use the last cached "current" index. (cherry picked from commit 5754e623261e992d965fccb25e3fcd81a59fa63c) --- .../inc/RooStats/HistFactory/ParamHistFunc.h | 2 +- roofit/histfactory/src/RooBarlowBeestonLL.cxx | 2 +- roofit/roofit/src/RooHistConstraint.cxx | 6 +-- roofit/roofit/src/RooIntegralMorph.cxx | 14 ++++--- roofit/roofit/src/RooLagrangianMorphFunc.cxx | 11 +++-- roofit/roofit/src/RooParamHistFunc.cxx | 9 ++-- roofit/roofitcore/src/RooAbsPdf.cxx | 22 +++++----- roofit/roofitcore/src/RooAbsReal.cxx | 2 +- roofit/roofitcore/src/RooBinnedGenContext.cxx | 19 +++++---- roofit/roofitcore/src/RooChi2Var.cxx | 4 +- roofit/roofitcore/src/RooDataHist.cxx | 16 ++++---- roofit/roofitcore/src/RooFFTConvPdf.cxx | 5 ++- roofit/roofitcore/src/RooHistFunc.cxx | 7 +--- roofit/roofitcore/src/RooHistPdf.cxx | 7 +--- roofit/roofitcore/test/testRooDataHist.cxx | 17 ++++---- roofit/roofitcore/test/testRooDataSet.cxx | 5 +-- roofit/roostats/src/MCMCInterval.cxx | 41 +++++++------------ roofit/xroofit/src/xRooNode.cxx | 6 +-- 18 files changed, 90 insertions(+), 105 deletions(-) diff --git a/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h b/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h index 07dbbfd48a1d9..7efad55d99417 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h @@ -48,7 +48,7 @@ class ParamHistFunc : public RooAbsReal { RooDataHist const& dataHist() const { return _dataSet; } - double binVolume() const { return _dataSet.binVolume(); } + double binVolume(std::size_t i) const { return _dataSet.binVolume(i); } bool forceAnalyticalInt(const RooAbsArg&) const override { return true ; } diff --git a/roofit/histfactory/src/RooBarlowBeestonLL.cxx b/roofit/histfactory/src/RooBarlowBeestonLL.cxx index d5bdb47cb6e4a..94251d16ef722 100644 --- a/roofit/histfactory/src/RooBarlowBeestonLL.cxx +++ b/roofit/histfactory/src/RooBarlowBeestonLL.cxx @@ -192,7 +192,7 @@ void RooStats::HistFactory::RooBarlowBeestonLL::initializeBarlowCache() { cache.bin_center = bin_center; cache.observables = obsSet; - cache.binVolume = param_func->binVolume(); + cache.binVolume = param_func->binVolume(bin_index); // Delete this part, simply a comment RooArgList obs_list( *(cache.bin_center) ); diff --git a/roofit/roofit/src/RooHistConstraint.cxx b/roofit/roofit/src/RooHistConstraint.cxx index fcff8e9307f4e..80cfc3f7113d9 100644 --- a/roofit/roofit/src/RooHistConstraint.cxx +++ b/roofit/roofit/src/RooHistConstraint.cxx @@ -60,11 +60,11 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, // Now populate nominal with parameters for (int i=0 ; i_dh.numEntries() ; i++) { - phf->_dh.get(i) ; - if (phf->_dh.weight()_dh.weight() != 0.) { + const double wi = phf->_dh.weight(i) ; + if (wi(vname,vname,0,1.E30); - var->setVal(phf->_dh.weight()) ; + var->setVal(wi) ; var->setConstant(true); auto gamma = static_cast(phf->_p.at(i)); diff --git a/roofit/roofit/src/RooIntegralMorph.cxx b/roofit/roofit/src/RooIntegralMorph.cxx index 9e2e1cc6cfdcb..1bc8519f7adee 100644 --- a/roofit/roofit/src/RooIntegralMorph.cxx +++ b/roofit/roofit/src/RooIntegralMorph.cxx @@ -392,8 +392,8 @@ void RooIntegralMorph::MorphCacheElem::calculate(TIterator* dIter) // Zero output histogram below lowest calculable X value for (int i=0; i<_yatXmin ; i++) { dIter->Next() ; - //_hist->get(i) ; - hist()->set(0) ; + const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true); + hist()->set(binIdx, 0, -1) ; } double x1 = _x->getMin("cache"); @@ -420,14 +420,16 @@ void RooIntegralMorph::MorphCacheElem::calculate(TIterator* dIter) double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ; dIter->Next() ; - //_hist->get(i) ; - hist()->set(fbarX) ; + { + const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true); + hist()->set(binIdx, fbarX, -1) ; + } } // Zero output histogram above highest calculable X value for (int i=_yatXmax+1 ; iNext() ; - //_hist->get(i) ; - hist()->set(0) ; + const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true); + hist()->set(binIdx, 0, -1) ; } pdf()->setUnitNorm(true) ; diff --git a/roofit/roofit/src/RooLagrangianMorphFunc.cxx b/roofit/roofit/src/RooLagrangianMorphFunc.cxx index c380ccb006221..94473baec0a1f 100644 --- a/roofit/roofit/src/RooLagrangianMorphFunc.cxx +++ b/roofit/roofit/src/RooLagrangianMorphFunc.cxx @@ -2572,12 +2572,12 @@ TH1 *RooLagrangianMorphFunc::createTH1(const std::string &name, bool correlateEr continue; } const RooDataHist &dhist = hf->dataHist(); - dhist.get(i); RooAbsReal *formula = dynamic_cast(prod->components().find(Form("w_%s", prod->GetName()))); double weight = formula->getVal(); - unc2 += dhist.weightSquared() * weight * weight; - unc += sqrt(dhist.weightSquared()) * weight; - val += dhist.weight() * weight; + const double w2 = dhist.weightSquared(i); + unc2 += w2 * weight * weight; + unc += sqrt(w2) * weight; + val += dhist.weight(i) * weight; } hist->SetBinContent(i + 1, val); hist->SetBinError(i + 1, correlateErrors ? unc : sqrt(unc2)); @@ -2827,8 +2827,7 @@ double RooLagrangianMorphFunc::expectedUncertainty() const if (hf) { const RooDataHist &hist = hf->dataHist(); for (Int_t j = 0; j < observable->getBins(); ++j) { - hist.get(j); - newunc2 += hist.weightSquared(); + newunc2 += hist.weightSquared(j); } } else if (rv) { newunc2 = pow(rv->getError(), 2); diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx index 2890c083c2c45..3f37cd548ece1 100644 --- a/roofit/roofit/src/RooParamHistFunc.cxx +++ b/roofit/roofit/src/RooParamHistFunc.cxx @@ -44,11 +44,11 @@ RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataH // Now populate p with parameters RooArgSet allVars; for (Int_t i = 0; i < _dh.numEntries(); i++) { - _dh.get(i); + const double wi = _dh.weight(i); const char *vname = Form("%s_gamma_bin_%i", GetName(), i); RooRealVar *var = new RooRealVar(vname, vname, 0, 1000); - var->setVal(_relParam ? 1 : _dh.weight()); - var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight())); + var->setVal(_relParam ? 1 : wi); + var->setError(_relParam ? 1 / sqrt(wi) : sqrt(wi)); var->setConstant(true); allVars.add(*var); _p.add(*var); @@ -94,8 +94,7 @@ void RooParamHistFunc::setActual(Int_t ibin, double newVal) double RooParamHistFunc::getNominal(Int_t ibin) const { - _dh.get(ibin) ; - return _dh.weight() ; + return _dh.weight(ibin) ; } //////////////////////////////////////////////////////////////////////////////// diff --git a/roofit/roofitcore/src/RooAbsPdf.cxx b/roofit/roofitcore/src/RooAbsPdf.cxx index 7aad2f55d3b9f..296649b6aae4e 100644 --- a/roofit/roofitcore/src/RooAbsPdf.cxx +++ b/roofit/roofitcore/src/RooAbsPdf.cxx @@ -1703,26 +1703,27 @@ RooFit::OwningPtr RooAbsPdf::generateBinned(const RooArgSet &whatVa Int_t histOutSum(0) ; for (int i=0 ; inumEntries() ; i++) { hist->get(i) ; + const double wi = hist->weight(i) ; if (expectedData) { // Expected data, multiply p.d.f by nEvents - double w=hist->weight()*nEvents ; + double w=wi*nEvents ; hist->set(i, w, sqrt(w)); } else if (extended) { // Extended mode, set contents to Poisson(pdf*nEvents) - double w = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ; - hist->set(w,sqrt(w)) ; + double w = RooRandom::randomGenerator()->Poisson(wi*nEvents) ; + hist->set(i, w, sqrt(w)) ; } else { // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill // histogram yet. - if (hist->weight()>histMax) { - histMax = hist->weight() ; + if (wi>histMax) { + histMax = wi ; } - histOut[i] = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ; + histOut[i] = RooRandom::randomGenerator()->Poisson(wi*nEvents) ; histOutSum += histOut[i] ; } } @@ -1745,7 +1746,7 @@ RooFit::OwningPtr RooAbsPdf::generateBinned(const RooArgSet &whatVa hist->get(ibinRand) ; double ranY = RooRandom::randomGenerator()->Uniform(histMax) ; - if (ranYweight()) { + if (ranYweight(ibinRand)) { if (wgt==1) { histOut[ibinRand]++ ; } else { @@ -1768,8 +1769,7 @@ RooFit::OwningPtr RooAbsPdf::generateBinned(const RooArgSet &whatVa // Transfer working array to histogram for (int i=0 ; inumEntries() ; i++) { - hist->get(i) ; - hist->set(histOut[i],sqrt(1.0*histOut[i])) ; + hist->set(i, histOut[i], sqrt(1.0*histOut[i])) ; } } else if (expectedData) { @@ -1779,8 +1779,8 @@ RooFit::OwningPtr RooAbsPdf::generateBinned(const RooArgSet &whatVa // bin average and bin integral in sampling bins double corr = nEvents/hist->sumEntries() ; for (int i=0 ; inumEntries() ; i++) { - hist->get(i) ; - hist->set(hist->weight()*corr,sqrt(hist->weight()*corr)) ; + const double wnew = hist->weight(i)*corr ; + hist->set(i, wnew, sqrt(wnew)) ; } } diff --git a/roofit/roofitcore/src/RooAbsReal.cxx b/roofit/roofitcore/src/RooAbsReal.cxx index b6e5a1f901749..da039f60240cd 100644 --- a/roofit/roofitcore/src/RooAbsReal.cxx +++ b/roofit/roofitcore/src/RooAbsReal.cxx @@ -1151,7 +1151,7 @@ RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSe const RooArgSet* obs = hist->get(i) ; double binVal = theClone->getVal(normSet?normSet:obs)*scaleFactor ; if (correctForBinSize) { - binVal*= hist->binVolume() ; + binVal*= hist->binVolume(i) ; } hist->set(i, binVal, 0.); } diff --git a/roofit/roofitcore/src/RooBinnedGenContext.cxx b/roofit/roofitcore/src/RooBinnedGenContext.cxx index dbcc5652e1bf3..e844a59dc6148 100644 --- a/roofit/roofitcore/src/RooBinnedGenContext.cxx +++ b/roofit/roofitcore/src/RooBinnedGenContext.cxx @@ -141,27 +141,28 @@ RooDataSet *RooBinnedGenContext::generate(double nEvt, bool /*skipInit*/, bool e double histMax(-1) ; Int_t histOutSum(0) ; for (int i=0 ; i<_hist->numEntries() ; i++) { - _hist->get(i) ; + const RooArgSet* coords = _hist->get(i) ; + const double wi = _hist->weight(i) ; if (_expectedData) { // Expected data, multiply p.d.f by nEvents - double w=_hist->weight()*nEvents ; - wudata->add(*_hist->get(),w) ; + double w=wi*nEvents ; + wudata->add(*coords,w) ; } else if (extended) { // Extended mode, set contents to Poisson(pdf*nEvents) - double w = RooRandom::randomGenerator()->Poisson(_hist->weight()*nEvents) ; - wudata->add(*_hist->get(),w) ; + double w = RooRandom::randomGenerator()->Poisson(wi*nEvents) ; + wudata->add(*coords,w) ; } else { // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill // histogram yet. - if (_hist->weight()>histMax) { - histMax = _hist->weight() ; + if (wi>histMax) { + histMax = wi ; } - histOut[i] = RooRandom::randomGenerator()->Poisson(_hist->weight()*nEvents) ; + histOut[i] = RooRandom::randomGenerator()->Poisson(wi*nEvents) ; histOutSum += histOut[i] ; } } @@ -182,7 +183,7 @@ RooDataSet *RooBinnedGenContext::generate(double nEvt, bool /*skipInit*/, bool e _hist->get(ibinRand) ; double ranY = RooRandom::randomGenerator()->Uniform(histMax) ; - if (ranY<_hist->weight()) { + if (ranY<_hist->weight(ibinRand)) { if (wgt==1) { histOut[ibinRand]++ ; } else { diff --git a/roofit/roofitcore/src/RooChi2Var.cxx b/roofit/roofitcore/src/RooChi2Var.cxx index b6df3e207403b..89020e32ece56 100644 --- a/roofit/roofitcore/src/RooChi2Var.cxx +++ b/roofit/roofitcore/src/RooChi2Var.cxx @@ -107,9 +107,9 @@ double RooChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEve } if (!doSelect) continue ; - const double nData = hdata->weight() ; + const double nData = hdata->weight(i) ; - const double nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume() ; + const double nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume(i) ; const double eExt = nPdf-nData ; diff --git a/roofit/roofitcore/src/RooDataHist.cxx b/roofit/roofitcore/src/RooDataHist.cxx index 5d6be15717fd8..c54b845d9fd88 100644 --- a/roofit/roofitcore/src/RooDataHist.cxx +++ b/roofit/roofitcore/src/RooDataHist.cxx @@ -660,7 +660,7 @@ void RooDataHist::importDHistSet(const RooArgList & /*vars*/, RooCategory &index // Transfer contents for (Int_t i = 0; i < dhist->numEntries(); i++) { _vars.assign(*dhist->get(i)); - add(_vars, dhist->weight() * initWgt, pow(dhist->weightError(SumW2), 2)); + add(_vars, dhist->weight(i) * initWgt, pow(dhist->weightError(SumW2), 2)); } } } @@ -955,7 +955,7 @@ std::unique_ptr RooDataHist::reduceEng(const RooArgSet& varSubset, c if (!cloneVar || cloneVar->getVal()) { weightError(lo,hi,SumW2) ; - rdh->add(*row,weight(),lo*lo) ; + rdh->add(*row,weight(i),lo*lo) ; } } @@ -1572,7 +1572,7 @@ void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Expected not allowed here"); break ; - case Poisson: + case Poisson: { if (_errLo && _errLo[_curIndex] >= 0.0) { // Weight is preset or precalculated lo = _errLo[_curIndex]; @@ -1586,12 +1586,14 @@ void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const // Calculate poisson errors double ym; double yp; - RooHistError::instance().getPoissonInterval(Int_t(weight()+0.5),ym,yp,1) ; - _errLo[_curIndex] = weight()-ym; - _errHi[_curIndex] = yp-weight(); + const double w = weight(_curIndex); + RooHistError::instance().getPoissonInterval(Int_t(w+0.5),ym,yp,1) ; + _errLo[_curIndex] = w-ym; + _errHi[_curIndex] = yp-w; lo = _errLo[_curIndex]; hi = _errHi[_curIndex]; return ; + } case SumW2: lo = std::sqrt(weightSquared(_curIndex)); @@ -2348,7 +2350,7 @@ void RooDataHist::printContents(std::ostream& os) const double lo, hi; weightError(lo, hi, RooAbsData::SumW2); - os << ", weight=" << weight() << " +/- [" << lo << "," << hi << "]" + os << ", weight=" << weight(i) << " +/- [" << lo << "," << hi << "]" << std::endl; } } diff --git a/roofit/roofitcore/src/RooFFTConvPdf.cxx b/roofit/roofitcore/src/RooFFTConvPdf.cxx index 15bfde3f3883a..790ea9aef1f5b 100644 --- a/roofit/roofitcore/src/RooFFTConvPdf.cxx +++ b/roofit/roofitcore/src/RooFFTConvPdf.cxx @@ -657,10 +657,11 @@ void RooFFTConvPdf::fillCacheSlice(FFTCacheElem& aux, const RooArgSet& slicePos) while (j>=N2) j-= N2 ; iter->Next() ; + const std::size_t binIdx = cacheHist.getIndex(*cacheHist.get(), /*fast=*/true); #ifndef ROOFIT_MATH_FFTW3 - cacheHist.set(output[j]); + cacheHist.set(binIdx, output[j], -1.); #else - cacheHist.set(aux.fftc2r->GetPointReal(j)); + cacheHist.set(binIdx, aux.fftc2r->GetPointReal(j), -1.); #endif } } diff --git a/roofit/roofitcore/src/RooHistFunc.cxx b/roofit/roofitcore/src/RooHistFunc.cxx index 54f9bf9e472cd..043df4cd492e4 100644 --- a/roofit/roofitcore/src/RooHistFunc.cxx +++ b/roofit/roofitcore/src/RooHistFunc.cxx @@ -247,8 +247,7 @@ double RooHistFunc::maxVal(Int_t code) const double max(-1) ; for (Int_t i=0 ; i<_dataHist->numEntries() ; i++) { - _dataHist->get(i) ; - double wgt = _dataHist->weight() ; + double wgt = _dataHist->weight(i) ; if (wgt>max) max=wgt ; } @@ -487,9 +486,7 @@ bool RooHistFunc::areIdentical(const RooDataHist& dh1, const RooDataHist& dh2) if (std::abs(dh1.sumEntries()-dh2.sumEntries())>1e-8) return false ; if (dh1.numEntries() != dh2.numEntries()) return false ; for (int i=0 ; i < dh1.numEntries() ; i++) { - dh1.get(i) ; - dh2.get(i) ; - if (std::abs(dh1.weight()-dh2.weight())>1e-8) return false ; + if (std::abs(dh1.weight(i)-dh2.weight(i))>1e-8) return false ; } using RooHelpers::getColonSeparatedNameString; if (getColonSeparatedNameString(*dh1.get()) != getColonSeparatedNameString(*dh2.get())) return false ; diff --git a/roofit/roofitcore/src/RooHistPdf.cxx b/roofit/roofitcore/src/RooHistPdf.cxx index dd3555827e854..4879f782596b3 100644 --- a/roofit/roofitcore/src/RooHistPdf.cxx +++ b/roofit/roofitcore/src/RooHistPdf.cxx @@ -530,8 +530,7 @@ double RooHistPdf::maxVal(Int_t code) const double max(-1) ; for (Int_t i=0 ; i<_dataHist->numEntries() ; i++) { - _dataHist->get(i) ; - double wgt = _dataHist->weight() ; + double wgt = _dataHist->weight(i) ; if (wgt>max) max=wgt ; } @@ -548,9 +547,7 @@ bool RooHistPdf::areIdentical(const RooDataHist& dh1, const RooDataHist& dh2) if (std::abs(dh1.sumEntries()-dh2.sumEntries())>1e-8) return false ; if (dh1.numEntries() != dh2.numEntries()) return false ; for (int i=0 ; i < dh1.numEntries() ; i++) { - dh1.get(i) ; - dh2.get(i) ; - if (std::abs(dh1.weight()-dh2.weight())>1e-8) return false ; + if (std::abs(dh1.weight(i)-dh2.weight(i))>1e-8) return false ; } return true ; } diff --git a/roofit/roofitcore/test/testRooDataHist.cxx b/roofit/roofitcore/test/testRooDataHist.cxx index e8d6084e75541..d2b78a81f94f3 100644 --- a/roofit/roofitcore/test/testRooDataHist.cxx +++ b/roofit/roofitcore/test/testRooDataHist.cxx @@ -96,7 +96,7 @@ TEST(RooDataHist, UnWeightedEntries) EXPECT_EQ(static_cast(coordsAtZero.find(x))->getVal(), static_cast(coordsAtPoint9.find(x))->getVal()); - const double weight = dataHist.weight(); + const double weight = dataHist.weight(dataHist.getIndex(coordinates)); EXPECT_EQ(weight, targetBinContent); EXPECT_NEAR(dataHist.weightError(RooAbsData::Poisson), sqrt(targetBinContent), @@ -124,7 +124,7 @@ TEST(RooDataHist, UnWeightedEntries) RooArgSet coordsAt10; dataHist.get(10)->snapshot(coordsAt10); - const double weightBin10 = dataHist.weight(); + const double weightBin10 = dataHist.weight(10); EXPECT_NEAR(static_cast(coordsAt10.find(x))->getVal(), 0.5, 1.E-1); EXPECT_EQ(weight, weightBin10); @@ -149,7 +149,7 @@ TEST(RooDataHist, WeightedEntries) x.setVal(0.); dataHist.get(coordinates); - const double weight = dataHist.weight(); + const double weight = dataHist.weight(dataHist.getIndex(coordinates)); ASSERT_EQ(weight, targetBinContent); const double targetError = sqrt(10 * 4.); @@ -178,7 +178,7 @@ TEST(RooDataHist, WeightedEntries) RooArgSet coordsAt10; dataHist.get(10)->snapshot(coordsAt10); - const double weightBin10 = dataHist.weight(); + const double weightBin10 = dataHist.weight(10); EXPECT_NEAR(static_cast(coordsAt10.find(x))->getVal(), 0.5, 1.E-1); EXPECT_EQ(weight, weightBin10); @@ -217,8 +217,7 @@ TEST_P(RooDataHistIO, ReadLegacy) EXPECT_EQ(dataHist.sumEntries(), 20 * targetBinContent); static_cast(legacyVals.find("x"))->setVal(0.); - dataHist.get(legacyVals); // trigger side effect for weight below. - const double weight = dataHist.weight(); + const double weight = dataHist.weight(dataHist.getIndex(legacyVals)); ASSERT_EQ(weight, targetBinContent); const double targetError = sqrt(10 * 4.); @@ -247,7 +246,7 @@ TEST_P(RooDataHistIO, ReadLegacy) RooArgSet coordsAt10; dataHist.get(10)->snapshot(coordsAt10); - const double weightBin10 = dataHist.weight(); + const double weightBin10 = dataHist.weight(10); EXPECT_NEAR(static_cast(coordsAt10.find("x"))->getVal(), 0.5, 1.E-1); EXPECT_EQ(weight, weightBin10); @@ -804,11 +803,11 @@ TEST(RooDataHist, SplitDataHistWithSumW2) data1.get(0); data2.get(0); - EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared()); + EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared(0)); EXPECT_DOUBLE_EQ(data2.weightError(), data1.weightError()); data1.get(1); data2.get(1); - EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared()); + EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared(1)); EXPECT_DOUBLE_EQ(data2.weightError(), data1.weightError()); } diff --git a/roofit/roofitcore/test/testRooDataSet.cxx b/roofit/roofitcore/test/testRooDataSet.cxx index c3d72473f4deb..7e1e45443feaa 100644 --- a/roofit/roofitcore/test/testRooDataSet.cxx +++ b/roofit/roofitcore/test/testRooDataSet.cxx @@ -374,10 +374,9 @@ TEST(RooDataSet, ImportDataHist) RooDataSet ds{"ds", "ds", x, RooFit::Import(dh)}; for (int i = 0; i < x.numBins(); ++i) { - dh.get(i); ds.get(i); - EXPECT_FLOAT_EQ(ds.weight(), dh.weight()) << "weight() is off in bin " << i; - EXPECT_FLOAT_EQ(ds.weightSquared(), dh.weightSquared()) << "weightSquared() is off in bin " << i; + EXPECT_FLOAT_EQ(ds.weight(), dh.weight(i)) << "weight() is off in bin " << i; + EXPECT_FLOAT_EQ(ds.weightSquared(), dh.weightSquared(i)) << "weightSquared() is off in bin " << i; } } diff --git a/roofit/roostats/src/MCMCInterval.cxx b/roofit/roostats/src/MCMCInterval.cxx index acfe5261503a9..3f246ccd79912 100644 --- a/roofit/roostats/src/MCMCInterval.cxx +++ b/roofit/roostats/src/MCMCInterval.cxx @@ -111,10 +111,8 @@ MCMCInterval::~MCMCInterval() = default; struct CompareDataHistBins { CompareDataHistBins(RooDataHist* hist) : fDataHist(hist) {} bool operator() (Int_t bin1 , Int_t bin2) { - fDataHist->get(bin1); - double n1 = fDataHist->weight(); - fDataHist->get(bin2); - double n2 = fDataHist->weight(); + double n1 = fDataHist->weight(bin1); + double n2 = fDataHist->weight(bin2); return (n1 < n2); } RooDataHist* fDataHist; @@ -185,8 +183,7 @@ bool MCMCInterval::IsInInterval(const RooArgSet& point) const // >= cutoff Int_t bin; bin = fDataHist->getIndex(point); - fDataHist->get(bin); - return (fDataHist->weight() >= (double)fHistCutoff); + return (fDataHist->weight(bin) >= (double)fHistCutoff); } } } else if (fIntervalType == kTailFraction) { @@ -839,8 +836,7 @@ void MCMCInterval::DetermineByDataHist() double content; Int_t i; for (i = numBins - 1; i >= 0; i--) { - fDataHist->get(bins[i]); - content = fDataHist->weight(); + content = fDataHist->weight(bins[i]); if ((sum + content) / nEntries >= fConfidenceLevel) { fHistCutoff = content; if (fIsHistStrict) { @@ -858,8 +854,7 @@ void MCMCInterval::DetermineByDataHist() if (fIsHistStrict) { // keep going to find the sum for ( ; i >= 0; i--) { - fDataHist->get(bins[i]); - content = fDataHist->weight(); + content = fDataHist->weight(bins[i]); if (content == fHistCutoff) { sum += content; } else { @@ -869,8 +864,7 @@ void MCMCInterval::DetermineByDataHist() } else { // backtrack to find the cutoff and sum for ( ; i < numBins; i++) { - fDataHist->get(bins[i]); - content = fDataHist->weight(); + content = fDataHist->weight(bins[i]); if (content > fHistCutoff) { fHistCutoff = content; break; @@ -1070,9 +1064,8 @@ double MCMCInterval::LowerLimitByDataHist(RooRealVar& param) double lowerLimit = param.getMax(); double val; for (Int_t i = 0; i < numBins; i++) { - fDataHist->get(i); - if (fDataHist->weight() >= fHistCutoff) { - val = fDataHist->get()->getRealValue(param.GetName()); + if (fDataHist->weight(i) >= fHistCutoff) { + val = fDataHist->get(i)->getRealValue(param.GetName()); if (val < lowerLimit) lowerLimit = val; } @@ -1147,9 +1140,8 @@ double MCMCInterval::UpperLimitByDataHist(RooRealVar& param) double upperLimit = param.getMin(); double val; for (Int_t i = 0; i < numBins; i++) { - fDataHist->get(i); - if (fDataHist->weight() >= fHistCutoff) { - val = fDataHist->get()->getRealValue(param.GetName()); + if (fDataHist->weight(i) >= fHistCutoff) { + val = fDataHist->get(i)->getRealValue(param.GetName()); if (val > upperLimit) upperLimit = val; } @@ -1187,9 +1179,8 @@ double MCMCInterval::LowerLimitByKeys(RooRealVar& param) double lowerLimit = param.getMax(); double val; for (Int_t i = 0; i < numBins; i++) { - fKeysDataHist->get(i); - if (fKeysDataHist->weight() >= fKeysCutoff) { - val = fKeysDataHist->get()->getRealValue(param.GetName()); + if (fKeysDataHist->weight(i) >= fKeysCutoff) { + val = fKeysDataHist->get(i)->getRealValue(param.GetName()); if (val < lowerLimit) lowerLimit = val; } @@ -1227,9 +1218,8 @@ double MCMCInterval::UpperLimitByKeys(RooRealVar& param) double upperLimit = param.getMin(); double val; for (Int_t i = 0; i < numBins; i++) { - fKeysDataHist->get(i); - if (fKeysDataHist->weight() >= fKeysCutoff) { - val = fKeysDataHist->get()->getRealValue(param.GetName()); + if (fKeysDataHist->weight(i) >= fKeysCutoff) { + val = fKeysDataHist->get(i)->getRealValue(param.GetName()); if (val > upperLimit) upperLimit = val; } @@ -1264,8 +1254,7 @@ double MCMCInterval::GetKeysMax() double max = 0; double w; for (Int_t i = 0; i < numBins; i++) { - fKeysDataHist->get(i); - w = fKeysDataHist->weight(); + w = fKeysDataHist->weight(i); if (w > max) max = w; } diff --git a/roofit/xroofit/src/xRooNode.cxx b/roofit/xroofit/src/xRooNode.cxx index 2bdeebf94c663..ec8d5034a8aa8 100644 --- a/roofit/xroofit/src/xRooNode.cxx +++ b/roofit/xroofit/src/xRooNode.cxx @@ -2136,7 +2136,7 @@ xRooNode xRooNode::Add(const xRooNode &child, Option_t *opt) // need to divide by bin widths first for (int i = 0; i < _f->dataHist().numEntries(); i++) { auto bin_pars = _f->dataHist().get(i); - _f->dataHist().set(*bin_pars, _f->dataHist().weight() / _f->dataHist().binVolume(*bin_pars)); + _f->dataHist().set(*bin_pars, _f->dataHist().weight(i) / _f->dataHist().binVolume(*bin_pars)); } _f->setAttribute("autodensity", false); _f->setValueDirty(); @@ -3329,7 +3329,7 @@ xRooNode xRooNode::Multiply(const xRooNode &child, Option_t *opt) // need to divide by bin widths first for (int i = 0; i < _f->dataHist().numEntries(); i++) { auto bin_pars = _f->dataHist().get(i); - _f->dataHist().set(*bin_pars, _f->dataHist().weight() / _f->dataHist().binVolume(*bin_pars)); + _f->dataHist().set(*bin_pars, _f->dataHist().weight(i) / _f->dataHist().binVolume(*bin_pars)); } _f->setValueDirty(); @@ -4893,7 +4893,7 @@ bool xRooNode::SetBinError(int bin, double value) TString origName = (f->getStringAttribute("origName")) ? f->getStringAttribute("origName") : GetName(); rrv->setStringAttribute(Form("sumw2_%s", origName.Data()), TString::Format("%f", pow(value, 2))); auto bin_pars = f->dataHist().get(bin - 1); - auto _binContent = f->dataHist().weight(); + auto _binContent = f->dataHist().weight(bin - 1); if (f->getAttribute("density")) { _binContent *= f->dataHist().binVolume(*bin_pars); } From ffb4333df69925143545b68808862d1481431447 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 21 Apr 2026 13:19:51 +0200 Subject: [PATCH 2/4] [RF][RS] Unify `BayesianCalculator::GetInterval()` 1D and ND code paths For models with a single POI and no nuisance parameters, `BayesianCalculator::GetInterval()` used `ComputeIntervalUsingRooFit()`, which built the CDF via `RooAbsPdf::createCdf()` On a computation graph that contains the NLL. The `createCdf()` relies on cloning the integrand subtree and substituting the POI with a prime variable, but the POI dependency inside the NLL compiled for the new vectorizing CPU evaluation backend is not propagated through the clone. The CDF and therefore the calculated interval is then wrong. Users worked around this by calling `SetScanOfPosterior(N)` instead (see #17567). Fortunately, there was already a code path for the n-dimensional parameter case present that doesn't use `createCdf()`, but the ROOT Math integrator classes directly. This commit commit suggests to also do this for the 1D POI-only case to avoid the issue with using the compiled NLL in `createCdf()` and simplifying the code by only providing one code path for interval caluclation with integration. Add a gtest regression test that reproduces the issue scenario and checks the default-path upper limit matches the explicit scan reference. Closes #17567. (cherry picked from commit 1eb9fb60899bc7b3a6a957147bc24dd44583c261) --- .../inc/RooStats/BayesianCalculator.h | 2 - roofit/roostats/src/BayesianCalculator.cxx | 98 ++++++------------- roofit/roostats/test/CMakeLists.txt | 1 + .../roostats/test/testBayesianCalculator.cxx | 65 ++++++++++++ 4 files changed, 98 insertions(+), 68 deletions(-) create mode 100644 roofit/roostats/test/testBayesianCalculator.cxx diff --git a/roofit/roostats/inc/RooStats/BayesianCalculator.h b/roofit/roostats/inc/RooStats/BayesianCalculator.h index f0800cb866abf..5529a456d362a 100644 --- a/roofit/roostats/inc/RooStats/BayesianCalculator.h +++ b/roofit/roostats/inc/RooStats/BayesianCalculator.h @@ -151,8 +151,6 @@ namespace RooStats { void ComputeIntervalFromCdf(double c1, double c2) const; - void ComputeIntervalUsingRooFit(double c1, double c2) const; - void ComputeShortestInterval() const; private: diff --git a/roofit/roostats/src/BayesianCalculator.cxx b/roofit/roostats/src/BayesianCalculator.cxx index d836e39948c75..3d5d2e513231d 100644 --- a/roofit/roostats/src/BayesianCalculator.cxx +++ b/roofit/roostats/src/BayesianCalculator.cxx @@ -54,7 +54,6 @@ credible interval from the given function. #include "RooAbsReal.h" #include "RooRealVar.h" #include "RooArgSet.h" -#include "RooBrentRootFinder.h" #include "RooFormulaVar.h" #include "RooGenericPdf.h" #include "RooPlot.h" @@ -167,7 +166,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fFunctor(nll, bindParams, RooArgList() ), // functor fPriorFunc(nullptr), fLikelihood(fFunctor, nullptr, nllMinimum), // integral of exp(-nll) function - fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ), // integrator fXmin(bindParams.size()), // vector of parameters (min values) fXmax(bindParams.size()) { @@ -176,7 +174,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fLikelihood.SetPrior(fPriorFunc.get() ); } - fIntegrator.SetFunction(fLikelihood, bindParams.size() ); + // ROOT::Math::IntegratorMultiDim does not support dim = 1, so fall back to + // the 1D integrator when there are no nuisance parameters. + if (bindParams.size() == 1) { + fIntegratorOneDim = std::make_unique(ROOT::Math::IntegratorOneDim::GetType(integType)); + fIntegratorOneDim->SetFunction(fLikelihood); + } else { + fIntegrator = std::make_unique(ROOT::Math::IntegratorMultiDim::GetType(integType)); + fIntegrator->SetFunction(fLikelihood, bindParams.size()); + } ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. " << " nllMinimum is " << nllMinimum << std::endl; @@ -191,7 +197,10 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { << " in interval [ " << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl; } - fIntegrator.Options().Print(ooccoutD(nullptr,NumericIntegration)); + if (fIntegrator) + fIntegrator->Options().Print(ooccoutD(nullptr,NumericIntegration)); + else + fIntegratorOneDim->Options().Print(ooccoutD(nullptr,NumericIntegration)); // store max POI value because it will be changed when evaluating the function fMaxPOI = fXmax[0]; @@ -214,7 +223,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { //fPriorFunc(std::shared_ptr((RooFunctor*)0)), fPriorFunc(rhs.fPriorFunc), fLikelihood(fFunctor, fPriorFunc.get(), rhs.fLikelihood.fOffset), - fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ), // integrator fXmin( rhs.fXmin), fXmax( rhs.fXmax), fNorm( rhs.fNorm), @@ -226,7 +234,13 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fError(rhs.fError), fNormCdfValues(rhs.fNormCdfValues) { - fIntegrator.SetFunction(fLikelihood, fXmin.size() ); + if (rhs.fIntegratorOneDim) { + fIntegratorOneDim = std::make_unique(ROOT::Math::IntegratorOneDim::GetType(rhs.fIntegratorOneDim->Name().c_str())); + fIntegratorOneDim->SetFunction(fLikelihood); + } else { + fIntegrator = std::make_unique(ROOT::Math::IntegratorMultiDim::GetType(rhs.fIntegrator->Name().c_str())); + fIntegrator->SetFunction(fLikelihood, fXmin.size()); + } // need special treatment for the smart pointer // if (rhs.fPriorFunc.get() ) { // fPriorFunc = std::shared_ptr(new RooFunctor(*(rhs.fPriorFunc) ) ); @@ -277,8 +291,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fFunctor.binding().resetNumCall(); // reset number of calls for debug - double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]); - double error = fIntegrator.Error(); + double cdf; + double error; + if (fIntegratorOneDim) { + cdf = fIntegratorOneDim->Integral(fXmin[0], fXmax[0]); + error = fIntegratorOneDim->Error(); + } else { + cdf = fIntegrator->Integral(&fXmin[0], &fXmax[0]); + error = fIntegrator->Error(); + } double normcdf = cdf/fNorm; // normalize the cdf ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , " @@ -324,7 +345,8 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { mutable RooFunctor fFunctor; // functor binding nll mutable std::shared_ptr fPriorFunc; // functor binding the prior LikelihoodFunction fLikelihood; // likelihood function - mutable ROOT::Math::IntegratorMultiDim fIntegrator; // integrator (mutable because Integral() is not const + mutable std::unique_ptr fIntegrator; // multi-dim integrator (for >=2 dims) + mutable std::unique_ptr fIntegratorOneDim; // 1D integrator (for POI-only case) mutable std::vector fXmin; // min value of parameters (poi+nuis) - mutable std::vector fXmax; // max value of parameters (poi+nuis) - max poi changes so it is mutable double fNorm = 1.0; // normalization value (computed in constructor) @@ -1127,14 +1149,7 @@ SimpleInterval* BayesianCalculator::GetInterval() const } else { - // use integration method if there are nuisance parameters - if (!fNuisanceParameters.empty()) { - ComputeIntervalFromCdf(lowerCutOff, upperCutOff); - } - else { - // case of no nuisance - just use createCdf from roofit - ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff); - } + ComputeIntervalFromCdf(lowerCutOff, upperCutOff); // case cdf failed (scan then the posterior) if (!fValidInterval) { fNScanBins = 100; @@ -1187,55 +1202,6 @@ double BayesianCalculator::GetMode() const { //return fApproxPosterior->GetMaximumX(); } -//////////////////////////////////////////////////////////////////////////////// -/// internal function compute the interval using RooFit - -void BayesianCalculator::ComputeIntervalUsingRooFit(double lowerCutOff, double upperCutOff ) const { - - coutI(Eval) << "BayesianCalculator: Compute interval using RooFit: posteriorPdf + createCdf + RooBrentRootFinder " << std::endl; - - RooRealVar* poi = dynamic_cast(fPOI.first() ); - assert(poi); - - fValidInterval = false; - if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf(); - if (!fPosteriorPdf) return; - - std::unique_ptr cdf{fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf())}; - if (!cdf) return; - - std::unique_ptr cdf_bind{cdf->bindVars(fPOI,&fPOI)}; - if (!cdf_bind) return; - - RooBrentRootFinder brf(*cdf_bind); - brf.setTol(fBrfPrecision); // set the brf precision - - double tmpVal = poi->getVal(); // patch used because findRoot changes the value of poi - - bool ret = true; - if (lowerCutOff > 0) { - double y = lowerCutOff; - ret &= brf.findRoot(fLower,poi->getMin(),poi->getMax(),y); - } - else - fLower = poi->getMin(); - - if (upperCutOff < 1.0) { - double y=upperCutOff; - ret &= brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y); - } - else - fUpper = poi->getMax(); - if (!ret) { - coutE(Eval) << "BayesianCalculator::GetInterval " - << "Error returned from Root finder, estimated interval is not fully correct" << std::endl; - } else { - fValidInterval = true; - } - - poi->setVal(tmpVal); // patch: restore the original value of poi -} - //////////////////////////////////////////////////////////////////////////////// /// internal function compute the interval using Cdf integration diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt index 75eff086a88d6..7770dc075e93e 100644 --- a/roofit/roostats/test/CMakeLists.txt +++ b/roofit/roostats/test/CMakeLists.txt @@ -1,4 +1,5 @@ ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats) +ROOT_ADD_GTEST(testBayesianCalculator testBayesianCalculator.cxx LIBRARIES RooStats) ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx LIBRARIES RooStats COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testHypoTestInvResult_1.root) diff --git a/roofit/roostats/test/testBayesianCalculator.cxx b/roofit/roostats/test/testBayesianCalculator.cxx new file mode 100644 index 0000000000000..164ed8a1357ae --- /dev/null +++ b/roofit/roostats/test/testBayesianCalculator.cxx @@ -0,0 +1,65 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +/// Regression test for https://github.com/root-project/root/issues/17567 + +TEST(BayesianCalculator, UpperLimitNoNuisance_Issue17567) +{ + RooRandom::randomGenerator()->SetSeed(12345); + + // Background-only generation model. + RooRealVar mass("mass", "mass", 6500, 7500); + RooRealVar slopeGen("slopeGen", "slopeGen", -8e-4); + RooExponential expoGen("expoGen", "expoGen", mass, slopeGen); + std::unique_ptr data{expoGen.generate(mass, 1000)}; + + // Fit model: Gaussian signal on top of exponential background. + // All parameters except the POI ngauss are fixed so there are no nuisances. + RooRealVar mean("mean", "mean", 7024); + RooRealVar sigma("sigma", "sigma", 44); + RooGaussian gauss("gauss", "gauss", mass, mean, sigma); + RooRealVar slope("slope", "slope", -8e-4); + RooExponential expo("expo", "expo", mass, slope); + RooRealVar ngauss("ngauss", "ngauss", 10, 0, 50); + RooRealVar nexpo("nexpo", "nexpo", 1e3); + RooAddPdf pdftot("pdftot", "pdftot", RooArgList(gauss, expo), RooArgList(ngauss, nexpo)); + + RooUniform prior("prior", "prior", ngauss); + + // Run with the default interval method (numerical integration instead of + // brute-force scan). + RooStats::BayesianCalculator bc(*data, pdftot, RooArgSet(ngauss), prior); + bc.SetConfidenceLevel(0.95); + bc.SetLeftSideTailFraction(0.); + std::unique_ptr bcInterval{bc.GetInterval()}; + ASSERT_NE(bcInterval, nullptr); + const double limitDefault = bcInterval->UpperLimit(); + EXPECT_EQ(bcInterval->LowerLimit(), 0.0); + + // Reference: explicit scan with many bins gives the same answer up to + // numerical integration and scan discretization errors. + RooStats::BayesianCalculator bcScan(*data, pdftot, RooArgSet(ngauss), prior); + bcScan.SetConfidenceLevel(0.95); + bcScan.SetLeftSideTailFraction(0.); + bcScan.SetScanOfPosterior(200); + std::unique_ptr scanInterval{bcScan.GetInterval()}; + ASSERT_NE(scanInterval, nullptr); + const double limitScan = scanInterval->UpperLimit(); + + // Both limits should sit well inside the POI range (if the old bug were + // present the limit would be at or near poi->getMax() = 50, or pulled down + // by a broken CDF root) and must agree within 1%. + EXPECT_GT(limitDefault, 20.0); + EXPECT_LT(limitDefault, 35.0); + EXPECT_NEAR(limitDefault, limitScan, 0.01 * limitScan); +} From 8d4f0bcc01cf97adae5600f8d31b0dfd597cadb1 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 20 Feb 2026 22:42:29 +0100 Subject: [PATCH 3/4] [RF] Disable redundant dirty-flag propagation during minimization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When a likelihood is evaluated with the new `"cpu"` backend, the `RooFit::Evaluator` fully manages dependency tracking and re-evaluation of the computation graph. In this case, RooFit’s built-in dirty flag propagation in RooAbsArg becomes redundant and introduces significant overhead for large models. This patch disables regular dirty state propagation for all non-fundamental nodes in the Evaluator's computation graph by setting their OperMode to `RooAbsArg::ADirty`. Fundamental nodes (e.g. RooRealVar, RooCategory) are excluded because they are often shared with other computation graphs outside the Evaluator (usually the original pdf in the RooWorkspace). To set the OperMode of *all* RooAbsArgs to `ADirty` during minimization, while avoiding side effects outside the minimization scope, the dirty flag propagation for the fundamental nodes is only disabled temporarily in the RooMinimizer. This commit drastically speeds up fits with AD in particular (up to 2 x for large models), because with fast gradients, the dirty flag propagation that determines which part of the compute graph needs to be recomputed becomes the bottleneck. It was also redundant with a faster "dirty state" bookkeeping mechanism in the `RooFit::Evaluator` class itself. At this point, there is no performance regression anymore when disabling recursive dirty flag propagation for all evaluated nodes, so the old comment in the code about test 14 in stressRooFit being slow doesn't apply anymore. (cherry picked from commit fa977749042b741c97375053703284a89a2ac1da) --- roofit/roofitcore/inc/RooEvaluatorWrapper.h | 3 +- roofit/roofitcore/inc/RooFit/Evaluator.h | 2 + roofit/roofitcore/inc/RooMinimizer.h | 1 + roofit/roofitcore/src/RooEvaluatorWrapper.cxx | 5 ++ roofit/roofitcore/src/RooFit/Evaluator.cxx | 64 ++++++++++++++++--- roofit/roofitcore/src/RooMinimizer.cxx | 48 ++++++++++---- 6 files changed, 102 insertions(+), 21 deletions(-) diff --git a/roofit/roofitcore/inc/RooEvaluatorWrapper.h b/roofit/roofitcore/inc/RooEvaluatorWrapper.h index dea80a4179602..340631e0d282d 100644 --- a/roofit/roofitcore/inc/RooEvaluatorWrapper.h +++ b/roofit/roofitcore/inc/RooEvaluatorWrapper.h @@ -71,9 +71,10 @@ class RooEvaluatorWrapper final : public RooAbsReal { void generateHessian(); void setUseGeneratedFunctionCode(bool); - void writeDebugMacro(std::string const &) const; + std::stack> setOperModes(RooAbsArg::OperMode opMode); + protected: double evaluate() const override; diff --git a/roofit/roofitcore/inc/RooFit/Evaluator.h b/roofit/roofitcore/inc/RooFit/Evaluator.h index 281e4ab7213f3..18bf0e304fac2 100644 --- a/roofit/roofitcore/inc/RooFit/Evaluator.h +++ b/roofit/roofitcore/inc/RooFit/Evaluator.h @@ -45,6 +45,8 @@ class Evaluator { void setOffsetMode(RooFit::EvalContext::OffsetMode); + std::stack> setOperModes(RooAbsArg::OperMode opMode); + private: void processVariable(NodeInfo &nodeInfo); void processCategory(NodeInfo &nodeInfo); diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index fcf2dfe960a64..95ff438588f3e 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -234,6 +234,7 @@ class RooMinimizer : public TObject { void fillCorrMatrix(RooFitResult &fitRes); void updateErrors(); + RooAbsReal &_function; ROOT::Fit::FitConfig _config; ///< fitter configuration (options and parameter settings) std::unique_ptr _result; /// _minimizer; ///writeDebugMacro(filename); } +std::stack> RooEvaluatorWrapper::setOperModes(RooAbsArg::OperMode opMode) +{ + return _evaluator->setOperModes(opMode); +} + } // namespace RooFit::Experimental /// \endcond diff --git a/roofit/roofitcore/src/RooFit/Evaluator.cxx b/roofit/roofitcore/src/RooFit/Evaluator.cxx index 0ad61d2774a26..87a7f929f76f6 100644 --- a/roofit/roofitcore/src/RooFit/Evaluator.cxx +++ b/roofit/roofitcore/src/RooFit/Evaluator.cxx @@ -46,6 +46,7 @@ RooAbsPdf::fitTo() is called and gets destroyed when the fitting ends. #include #include #include +#include namespace RooFit { @@ -325,16 +326,16 @@ void Evaluator::updateOutputSizes() for (auto &info : _nodes) { info.outputSize = outputSizeMap.at(info.absArg); info.isDirty = true; - - // In principle we don't need dirty flag propagation because the driver - // takes care of deciding which node needs to be re-evaluated. However, - // disabling it also for scalar mode results in very long fitting times - // for specific models (test 14 in stressRooFit), which still needs to be - // understood. TODO. - if (!info.isScalar()) { + // We don't need dirty flag propagation because the evaluator takes care + // of deciding what needs to be re-evaluated. We can disable the regular + // dirty state propagation. However, fundamental variables like + // RooRealVars and RooCategories are usually shared with other + // computation graphs outside the evaluator, so they can't be mutated. + // See also the code of the RooMinimizer, which ensures that dirty state + // propagation is temporarily disabled during minimization to really + // eliminate any overhead from the dirty flag propagation. + if (!info.absArg->isFundamental()) { setOperMode(info.absArg, RooAbsArg::ADirty); - } else { - setOperMode(info.absArg, info.originalOperMode); } } @@ -632,6 +633,51 @@ void Evaluator::setOperMode(RooAbsArg *arg, RooAbsArg::OperMode opMode) } } +// Change the operation modes of all RooAbsArgs in the computation graph. +// The changes are reset when the returned RAII object goes out of scope. +// +// We also walk transitively through value clients of the nodes to cover any +// node that RooAbsReal::doEval (the fallback scalar implementation) might +// inadvertently propagate the ADirty mode to via its recursive restore: that +// helper sets servers temporarily to AClean and then calls +// setOperMode(oldOperMode) to restore, which recurses to value clients when +// oldOperMode is ADirty. If we did not protect those clients here, any node +// outside the computation graph that shares a fundamental (e.g. a parameter +// like a RooRealVar) would be left permanently in ADirty after the first +// minimization, dramatically slowing down later scalar evaluations (for +// example on pdfs held by the legacy test statistics' internal cache). +std::stack> Evaluator::setOperModes(RooAbsArg::OperMode opMode) +{ + std::stack> out{}; + std::unordered_set visited; + + std::vector queue; + queue.reserve(_nodes.size()); + for (auto &info : _nodes) { + queue.push_back(info.absArg); + } + + while (!queue.empty()) { + RooAbsArg *node = queue.back(); + queue.pop_back(); + if (!visited.insert(node).second) + continue; + + if (opMode != node->operMode()) { + out.emplace(std::make_unique(node, opMode)); + } + + // Only follow value-client links: that is exactly the propagation path + // used by RooAbsArg::setOperMode with mode==ADirty. + if (opMode == RooAbsArg::ADirty) { + for (auto *client : node->valueClients()) { + queue.push_back(client); + } + } + } + return out; +} + void Evaluator::print(std::ostream &os) { std::cout << "--- RooFit BatchMode evaluation ---\n"; diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 718a13880ac2b..679a7c98a0450 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -40,27 +40,30 @@ automatic PDF optimization. #include "RooMinimizer.h" #include "RooAbsMinimizerFcn.h" -#include "RooArgSet.h" -#include "RooArgList.h" #include "RooAbsReal.h" +#include "RooArgList.h" +#include "RooArgSet.h" +#include "RooCategory.h" #include "RooDataSet.h" -#include "RooRealVar.h" -#include "RooSentinel.h" +#include "RooEvaluatorWrapper.h" +#include "RooFit/TestStatistics/RooAbsL.h" +#include "RooFit/TestStatistics/RooRealL.h" +#include "RooFitResult.h" +#include "RooHelpers.h" +#include "RooMinimizerFcn.h" #include "RooMsgService.h" -#include "RooCategory.h" #include "RooMultiPdf.h" #include "RooPlot.h" -#include "RooHelpers.h" -#include "RooMinimizerFcn.h" -#include "RooFitResult.h" -#include "RooFit/TestStatistics/RooAbsL.h" -#include "RooFit/TestStatistics/RooRealL.h" +#include "RooRealVar.h" +#include "RooSentinel.h" #ifdef ROOFIT_MULTIPROCESS #include "TestStatistics/MinuitFcnGrad.h" #include "RooFit/MultiProcess/Config.h" #include "RooFit/MultiProcess/ProcessTimer.h" #endif +#include "RooFitImplHelpers.h" + #include #include #include @@ -120,6 +123,22 @@ void reorderCombinations(std::vector> &combos, const std::vecto } } +// The RooEvaluatorWrapper uses its own logic to decide what needs to be +// re-evaluated. We can therefore disable the regular dirty state propagation +// temporarily during minimization. However, some RooAbsArgs shared with other +// regular RooFit computation graphs outside the minimized likelihood, so we +// have to make sure that the operation mode is reset after the minimization. +// +// This should be called before running any routine via the _minimizer data +// member. The RAII object should only be destructed after the routine is done. +std::stack> setOperModesDirty(RooAbsReal &function) +{ + if (auto *wrapper = dynamic_cast(&function)) { + return wrapper->setOperModes(RooAbsArg::ADirty); + } + return {}; +} + } // namespace //////////////////////////////////////////////////////////////////////////////// @@ -135,7 +154,7 @@ void reorderCombinations(std::vector> &combos, const std::vecto /// value of the input function. /// Constructor that accepts all configuration in struct with RooAbsReal likelihood -RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _cfg(cfg) +RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _function{function}, _cfg(cfg) { initMinimizerFirstPart(); auto nll_real = dynamic_cast(&function); @@ -692,6 +711,7 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do n[4] = n5; n[5] = n6; + auto operModeRAII = setOperModesDirty(_function); for (int ic = 0; ic < 6; ic++) { if (n[ic] > 0) { @@ -906,6 +926,8 @@ bool RooMinimizer::fitFCN() // fit a user provided FCN function // create fit parameter settings + auto operModeRAII = setOperModesDirty(_function); + // Check number of parameters unsigned int npar = getNPar(); if (npar == 0) { @@ -1045,6 +1067,8 @@ bool RooMinimizer::calculateHessErrors() // compute the Hesse errors according to configuration // set in the parameters and append value in fit result + auto operModeRAII = setOperModesDirty(_function); + // update minimizer (recreate if not done or if name has changed if (!updateMinimizerOptions()) { coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; @@ -1079,6 +1103,8 @@ bool RooMinimizer::calculateMinosErrors() // (in DoMinimization) aftewr minimizing if the // FitConfig::MinosErrors() flag is set + auto operModeRAII = setOperModesDirty(_function); + // update minimizer (but cannot re-create in this case). Must use an existing one if (!updateMinimizerOptions(false)) { coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl; From 8837e4e9198af4a9f8eee9804922768b117521a1 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 23 Apr 2026 00:03:31 +0200 Subject: [PATCH 4/4] [RF] Refactor ChangeOperModeRAII to work with groups of RooAbsArgs Several places needed to record a set of operation-mode changes and restore them later as a group, so it's better to have the ChangeOperModeRAII act on groups of RooAbsArg to not have to create one RAII object per arg. (cherry picked from commit d42a27ca7b0972c08b471eca650f12028b7f3291) --- roofit/roofitcore/inc/RooEvaluatorWrapper.h | 4 +- roofit/roofitcore/inc/RooFit/Evaluator.h | 5 +- roofit/roofitcore/res/RooFitImplHelpers.h | 48 ++++++++++++++----- roofit/roofitcore/src/RooEvaluatorWrapper.cxx | 2 +- roofit/roofitcore/src/RooFit/Evaluator.cxx | 14 +++--- roofit/roofitcore/src/RooMinimizer.cxx | 2 +- roofit/roofitcore/src/RooRealIntegral.cxx | 4 +- 7 files changed, 52 insertions(+), 27 deletions(-) diff --git a/roofit/roofitcore/inc/RooEvaluatorWrapper.h b/roofit/roofitcore/inc/RooEvaluatorWrapper.h index 340631e0d282d..ce394fc629152 100644 --- a/roofit/roofitcore/inc/RooEvaluatorWrapper.h +++ b/roofit/roofitcore/inc/RooEvaluatorWrapper.h @@ -21,8 +21,10 @@ #include #include +#include #include +class ChangeOperModeRAII; class RooAbsArg; class RooAbsCategory; class RooAbsPdf; @@ -73,7 +75,7 @@ class RooEvaluatorWrapper final : public RooAbsReal { void setUseGeneratedFunctionCode(bool); void writeDebugMacro(std::string const &) const; - std::stack> setOperModes(RooAbsArg::OperMode opMode); + std::unique_ptr setOperModes(RooAbsArg::OperMode opMode); protected: double evaluate() const override; diff --git a/roofit/roofitcore/inc/RooFit/Evaluator.h b/roofit/roofitcore/inc/RooFit/Evaluator.h index 18bf0e304fac2..f50424c634abf 100644 --- a/roofit/roofitcore/inc/RooFit/Evaluator.h +++ b/roofit/roofitcore/inc/RooFit/Evaluator.h @@ -20,7 +20,6 @@ #include #include -#include class ChangeOperModeRAII; class RooAbsArg; @@ -45,7 +44,7 @@ class Evaluator { void setOffsetMode(RooFit::EvalContext::OffsetMode); - std::stack> setOperModes(RooAbsArg::OperMode opMode); + std::unique_ptr setOperModes(RooAbsArg::OperMode opMode); private: void processVariable(NodeInfo &nodeInfo); @@ -68,7 +67,7 @@ class Evaluator { RooFit::EvalContext _evalContextCUDA; std::vector _nodes; // the ordered computation graph std::unordered_map _nodesMap; // for quick lookup of nodes - std::stack> _changeOperModeRAIIs; + std::unique_ptr _operModeChanges; }; } // end namespace RooFit diff --git a/roofit/roofitcore/res/RooFitImplHelpers.h b/roofit/roofitcore/res/RooFitImplHelpers.h index 6d04134bc85ec..d14bc905cd21e 100644 --- a/roofit/roofitcore/res/RooFitImplHelpers.h +++ b/roofit/roofitcore/res/RooFitImplHelpers.h @@ -46,25 +46,51 @@ class DisableCachingRAII { bool _oldState; }; -/// Struct to temporarily change the operation mode of a RooAbsArg until it -/// goes out of scope. +/// Scope guard that temporarily changes the operation mode of one or more +/// RooAbsArg instances. Each call to change() records the arg's current +/// operMode before flipping it to the requested mode (non-recursively, i.e. +/// value clients are not touched). Destruction (or an explicit clear()) +/// restores every recorded mode in LIFO order. +/// +/// The class is movable but not copyable, so it can be returned from +/// functions that build up a batch of changes to hand to the caller. class ChangeOperModeRAII { public: - ChangeOperModeRAII(RooAbsArg *arg, RooAbsArg::OperMode opMode) : _arg{arg}, _oldOpMode(arg->operMode()) + ChangeOperModeRAII() = default; + + /// Convenience ctor: behaves like a scope guard for a single arg. + ChangeOperModeRAII(RooAbsArg *arg, RooAbsArg::OperMode opMode) { change(arg, opMode); } + + ~ChangeOperModeRAII() { clear(); } + + ChangeOperModeRAII(ChangeOperModeRAII &&) = default; + ChangeOperModeRAII &operator=(ChangeOperModeRAII &&) = default; + ChangeOperModeRAII(ChangeOperModeRAII const &) = delete; + ChangeOperModeRAII &operator=(ChangeOperModeRAII const &) = delete; + + /// Record arg's current operMode and flip it to opMode. If the current + /// mode already equals opMode, this is a no-op (nothing to restore). + void change(RooAbsArg *arg, RooAbsArg::OperMode opMode) { - arg->setOperMode(opMode, /*recurse=*/false); + if (opMode == arg->operMode()) + return; + _entries.emplace_back(arg, arg->operMode()); + arg->setOperMode(opMode, /*recurseADirty=*/false); } - ChangeOperModeRAII(ChangeOperModeRAII const &other) = delete; - ChangeOperModeRAII &operator=(ChangeOperModeRAII const &other) = delete; - ChangeOperModeRAII(ChangeOperModeRAII &&other) = delete; - ChangeOperModeRAII &operator=(ChangeOperModeRAII &&other) = delete; + /// Restore every recorded change right away, emptying this guard. + void clear() + { + for (auto it = _entries.rbegin(); it != _entries.rend(); ++it) { + it->first->setOperMode(it->second, /*recurseADirty=*/false); + } + _entries.clear(); + } - ~ChangeOperModeRAII() { _arg->setOperMode(_oldOpMode, /*recurse=*/false); } + bool empty() const { return _entries.empty(); } private: - RooAbsArg *_arg = nullptr; - RooAbsArg::OperMode _oldOpMode; + std::vector> _entries; }; namespace RooHelpers { diff --git a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx index f3916c3e9e655..c0a2a66c0b1e5 100644 --- a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx +++ b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx @@ -741,7 +741,7 @@ void RooEvaluatorWrapper::writeDebugMacro(std::string const &filename) const return _funcWrapper->writeDebugMacro(filename); } -std::stack> RooEvaluatorWrapper::setOperModes(RooAbsArg::OperMode opMode) +std::unique_ptr RooEvaluatorWrapper::setOperModes(RooAbsArg::OperMode opMode) { return _evaluator->setOperModes(opMode); } diff --git a/roofit/roofitcore/src/RooFit/Evaluator.cxx b/roofit/roofitcore/src/RooFit/Evaluator.cxx index 87a7f929f76f6..b17fc94cbe598 100644 --- a/roofit/roofitcore/src/RooFit/Evaluator.cxx +++ b/roofit/roofitcore/src/RooFit/Evaluator.cxx @@ -628,9 +628,9 @@ void Evaluator::markGPUNodes() /// Evaluator gets deleted. void Evaluator::setOperMode(RooAbsArg *arg, RooAbsArg::OperMode opMode) { - if (opMode != arg->operMode()) { - _changeOperModeRAIIs.emplace(std::make_unique(arg, opMode)); - } + if (!_operModeChanges) + _operModeChanges = std::make_unique(); + _operModeChanges->change(arg, opMode); } // Change the operation modes of all RooAbsArgs in the computation graph. @@ -646,9 +646,9 @@ void Evaluator::setOperMode(RooAbsArg *arg, RooAbsArg::OperMode opMode) // like a RooRealVar) would be left permanently in ADirty after the first // minimization, dramatically slowing down later scalar evaluations (for // example on pdfs held by the legacy test statistics' internal cache). -std::stack> Evaluator::setOperModes(RooAbsArg::OperMode opMode) +std::unique_ptr Evaluator::setOperModes(RooAbsArg::OperMode opMode) { - std::stack> out{}; + auto out = std::make_unique(); std::unordered_set visited; std::vector queue; @@ -663,9 +663,7 @@ std::stack> Evaluator::setOperModes(RooAbsAr if (!visited.insert(node).second) continue; - if (opMode != node->operMode()) { - out.emplace(std::make_unique(node, opMode)); - } + out->change(node, opMode); // Only follow value-client links: that is exactly the propagation path // used by RooAbsArg::setOperMode with mode==ADirty. diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 679a7c98a0450..af81c50102902 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -131,7 +131,7 @@ void reorderCombinations(std::vector> &combos, const std::vecto // // This should be called before running any routine via the _minimizer data // member. The RAII object should only be destructed after the routine is done. -std::stack> setOperModesDirty(RooAbsReal &function) +std::unique_ptr setOperModesDirty(RooAbsReal &function) { if (auto *wrapper = dynamic_cast(&function)) { return wrapper->setOperModes(RooAbsArg::ADirty); diff --git a/roofit/roofitcore/src/RooRealIntegral.cxx b/roofit/roofitcore/src/RooRealIntegral.cxx index c580bb0b2c10e..a1aa4bc459e54 100644 --- a/roofit/roofitcore/src/RooRealIntegral.cxx +++ b/roofit/roofitcore/src/RooRealIntegral.cxx @@ -841,12 +841,12 @@ double RooRealIntegral::evaluate() const // from caching subgraph results (e.g. for nested numeric integrals). RooArgList serverList; _function->treeNodeServerList(&serverList, nullptr, true, true, false, true); - std::list operModeRAII; + ChangeOperModeRAII operModeRAII; for (auto *arg : serverList) { arg->syncCache(); if (arg->operMode() == RooAbsArg::AClean) { - operModeRAII.emplace_back(arg, RooAbsArg::Auto); + operModeRAII.change(arg, RooAbsArg::Auto); } }