From cfbff6a0e538d81484b050476fb3a92caa47b182 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sat, 1 Feb 2025 13:19:16 +0100 Subject: [PATCH 1/2] [RF] Correctly override `RooWrapperPdf::forceAnalyticalInt()` Just like with other wrapper classes like RooExtendPdf, we can safely use the analytical integration capabilities of the wrapped object, becuase we don't do any no-linear transformation. --- roofit/roofitcore/inc/RooWrapperPdf.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/roofit/roofitcore/inc/RooWrapperPdf.h b/roofit/roofitcore/inc/RooWrapperPdf.h index b1fb310f56744..aab16f468e8df 100644 --- a/roofit/roofitcore/inc/RooWrapperPdf.h +++ b/roofit/roofitcore/inc/RooWrapperPdf.h @@ -53,8 +53,11 @@ class RooWrapperPdf final : public RooAbsPdf { bool selfNormalized() const override { return _selfNormalized; } // Analytical Integration handling - bool forceAnalyticalInt(const RooAbsArg& dep) const override { - return _func->forceAnalyticalInt(dep); + bool forceAnalyticalInt(const RooAbsArg& /*dep*/) const override { + // Just like with other wrapper classes like RooExtendPdf, we can safely + // use the analytical integration capabilities of the wrapped object, + // becuase we don't do any no-linear transformation. + return true; } Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet, const char* rangeName=nullptr) const override { From f2da82cdb753ad9b4543db384c27ea97db6130ec Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sat, 1 Feb 2025 13:21:16 +0100 Subject: [PATCH 2/2] [RF] Support analytical integration of RooStepFunction Implement analytical integration capability for RooStepFunction, because this class is the easiest way to implement a bounded uniform distribution without adding serialized class members to the RooUniform (see #7880). Add a unit test that verifies that wrapping a `RooStepFunction` in a `RooWrapperPdf` doesn't result in any numerical integrals. --- roofit/roofit/inc/RooStepFunction.h | 5 ++ roofit/roofit/src/RooStepFunction.cxx | 64 ++++++++++++------- roofit/roofit/test/CMakeLists.txt | 1 + roofit/roofit/test/testRooStepFunc.cxx | 35 ++++++++++ .../roofitcore/inc/RooFit/Detail/MathFuncs.h | 14 ++++ 5 files changed, 97 insertions(+), 22 deletions(-) create mode 100644 roofit/roofit/test/testRooStepFunc.cxx diff --git a/roofit/roofit/inc/RooStepFunction.h b/roofit/roofit/inc/RooStepFunction.h index 879e660a88904..58714b689a159 100644 --- a/roofit/roofit/inc/RooStepFunction.h +++ b/roofit/roofit/inc/RooStepFunction.h @@ -38,6 +38,9 @@ class RooStepFunction : public RooAbsReal { std::list* plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const override ; + int getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=nullptr) const override; + double analyticalIntegral(int code, const char* rangeName=nullptr) const override; + protected: double evaluate() const override; @@ -48,6 +51,8 @@ class RooStepFunction : public RooAbsReal { RooListProxy _coefList ; RooListProxy _boundaryList ; bool _interpolate = false; + mutable std::vector _coefCache; //! + mutable std::vector _boundaryCache; //! ClassDefOverride(RooStepFunction,1) // Step Function }; diff --git a/roofit/roofit/src/RooStepFunction.cxx b/roofit/roofit/src/RooStepFunction.cxx index 6c9691010672c..fb9bec4ff658d 100644 --- a/roofit/roofit/src/RooStepFunction.cxx +++ b/roofit/roofit/src/RooStepFunction.cxx @@ -32,12 +32,26 @@ but a not-normalized function (RooAbsReal) #include #include -#include +#include #include -#include +#include ClassImp(RooStepFunction); +namespace { + +// Get the values for each element in a RooListProxy, using the correct normalization set. +std::span values(RooListProxy const &listProxy, std::vector &valueCache) +{ + valueCache.resize(listProxy.size()); + for (std::size_t i = 0; i < listProxy.size(); ++i) { + valueCache[i] = static_cast(listProxy[i]).getVal(listProxy.nset()); + } + return valueCache; +} + +} // namespace + //////////////////////////////////////////////////////////////////////////////// /// Constructor @@ -77,12 +91,8 @@ RooStepFunction::RooStepFunction(const RooStepFunction& other, const char* name) double RooStepFunction::evaluate() const { - std::vector b(_boundaryList.size()) ; - std::vector c(_coefList.size()+3) ; - Int_t nb(0) ; - for (auto * boundary : static_range_cast(_boundaryList)) { - b[nb++] = boundary->getVal() ; - } + std::span b = values(_boundaryList, _boundaryCache); + int nb = b.size(); // Return zero if outside any boundaries if ((_xb[nb-1])) return 0 ; @@ -90,25 +100,27 @@ double RooStepFunction::evaluate() const if (!_interpolate) { // No interpolation -- Return values bin-by-bin - for (Int_t i=0;ib[i]&&_x<=b[i+1]) { return (static_cast(_coefList.at(i)))->getVal() ; } } return 0 ; - } else { + } + + std::vector c(_coefList.size()+3) ; // Interpolation // Make array of (b[0],bin centers,b[last]) c[0] = b[0] ; c[nb] = b[nb-1] ; - for (Int_t i=0 ; i y(_coefList.size()+3) ; y[nc++] = 0 ; for(auto * coef : static_range_cast(_coefList)) { @@ -116,7 +128,7 @@ double RooStepFunction::evaluate() const } y[nc++] = 0 ; - for (Int_t i=0;ic[i]&&_x<=c[i+1]) { double xx[2] ; xx[0]=c[i] ; xx[1]=c[i+1] ; double yy[2] ; yy[0]=y[i] ; yy[1]=y[i+1] ; @@ -124,7 +136,6 @@ double RooStepFunction::evaluate() const } } return 0; - } } @@ -134,14 +145,23 @@ std::list *RooStepFunction::plotSamplingHint(RooAbsRealLValue &obs, doub return nullptr; } - // Retrieve position of all bin boundaries - std::vector boundaries; - boundaries.reserve(_boundaryList.size()); - for (auto *boundary : static_range_cast(_boundaryList)) { - boundaries.push_back(boundary->getVal()); - } - // Use the helper function from RooCurve to make sure to get sampling hints // that work with the RooFitPlotting. - return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi); + return RooCurve::plotSamplingHintForBinBoundaries(values(_boundaryList, _boundaryCache), xlo, xhi); +} + +int RooStepFunction::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const +{ + // We only support analytical integration if we integrate over x and there is + // no interpolation. + if (!_interpolate && matchArgs(allVars, analVars, _x)) + return 1; + return 0; +} + +double RooStepFunction::analyticalIntegral(int /*code*/, const char *rangeName) const +{ + return RooFit::Detail::MathFuncs::stepFunctionIntegral(_x.min(rangeName), _x.max(rangeName), _coefList.size(), + values(_boundaryList, _boundaryCache).data(), + values(_coefList, _coefCache).data()); } diff --git a/roofit/roofit/test/CMakeLists.txt b/roofit/roofit/test/CMakeLists.txt index 7169922134c61..a07a4adc1cb43 100644 --- a/roofit/roofit/test/CMakeLists.txt +++ b/roofit/roofit/test/CMakeLists.txt @@ -19,6 +19,7 @@ ROOT_ADD_GTEST(testRooKeysPdf testRooKeysPdf.cxx LIBRARIES Gpad RooFit) ROOT_ADD_GTEST(testRooLandau testRooLandau.cxx LIBRARIES RooFit) ROOT_ADD_GTEST(testRooParamHistFunc testRooParamHistFunc.cxx LIBRARIES Gpad RooFit) ROOT_ADD_GTEST(testRooPoisson testRooPoisson.cxx LIBRARIES RooFit) +ROOT_ADD_GTEST(testRooStepFunc testRooStepFunc.cxx LIBRARIES RooFit) if(mathmore) ROOT_EXECUTABLE(testRooFit testRooFit.cxx LIBRARIES RooFit MathMore) diff --git a/roofit/roofit/test/testRooStepFunc.cxx b/roofit/roofit/test/testRooStepFunc.cxx new file mode 100644 index 0000000000000..e24c0eb6c339e --- /dev/null +++ b/roofit/roofit/test/testRooStepFunc.cxx @@ -0,0 +1,35 @@ +// Tests for the RooStepFunc +// Authors: Jonas Rembser, CERN 2025 + +#include +#include +#include +#include + +#include + +// Test that RooStepFunction can be wrapped into a pdf object, where it is +// normalized without any analytical integrals. +TEST(RooStepFunc, InPdfWrapper) +{ + RooRealVar x{"x", "x", 5.0, 0.0, 10.0}; + RooRealVar a{"a", "a", 3.0, 0.0, 10.0}; + RooRealVar b{"b", "b", 7.0, 0.0, 10.0}; + + RooStepFunction stepFunc{"step", "", x, RooArgList{1.0}, {a, b}}; + RooWrapperPdf stepPdf{"pdf", "", stepFunc}; + + RooArgSet normSet{x}; + + RooArgSet intSet{x}; + RooArgSet numSet{}; + + RooRealIntegral integ{"integ", "", stepPdf, x}; + + // Make sure that the normalization integral for the pdf is analytical + EXPECT_EQ(integ.anaIntVars().size(), 1); + EXPECT_DOUBLE_EQ(integ.getVal(), 4.0); + + // Check that the pdf in correctly normalized + EXPECT_DOUBLE_EQ(stepPdf.getVal(normSet), 0.25); +} diff --git a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h index 22865c1f33fd9..07299e1ca5bb7 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h +++ b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h @@ -765,6 +765,20 @@ inline double multiVarGaussian(int n, const double *x, const double *mu, const d return std::exp(-0.5 * result); } +// Integral of a step function defined by `nBins` intervals, where the +// intervals have values `coefs` and the boundary on the interval `iBin` is +// given by `[boundaries[i], boundaries[i+1])`. +inline double stepFunctionIntegral(double xmin, double xmax, std::size_t nBins, double const *boundaries, double const *coefs) +{ + double out = 0.0; + for (std::size_t i = 0; i < nBins; ++i) { + double a = boundaries[i]; + double b = boundaries[i + 1]; + out += coefs[i] * std::max(0.0, std::min(b, xmax) - std::max(a, xmin)); + } + return out; +} + } // namespace MathFuncs } // namespace Detail } // namespace RooFit