Skip to content

[RF] Generalize boundary finding for RooHistPdf and RooHistFunc #21969

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
guitargeek wants to merge 1 commit into
base: master
Choose a base branch
from
+197 −123
Draft

[RF] Generalize boundary finding for RooHistPdf and RooHistFunc #21969

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
5 changes: 5 additions & 0 deletions roofit/roofitcore/inc/RooHistPdf.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include "RooDataHist.h"

#include <list>
#include <vector>

class RooRealVar;
class RooAbsReal;
Expand Down Expand Up @@ -145,6 +146,10 @@ class RooHistPdf : public RooAbsPdf {
double xlo,
double xhi);

static std::vector<double> histogramBoundariesInPlotObs(RooDataHist const &dataHist, RooArgSet const &pdfObsList,
RooArgSet const &histObsList, RooAbsRealLValue &obs,
double xlo, double xhi);

inline void initializeOwnedDataHist(std::unique_ptr<RooDataHist> &&dataHist)
{
_ownedDataHist = std::move(dataHist);
Expand Down
95 changes: 15 additions & 80 deletions roofit/roofitcore/src/RooHistFunc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -326,11 +326,11 @@ std::list<double>* RooHistFunc::plotSamplingHint(RooAbsRealLValue& obs, double x
return RooHistPdf::plotSamplingHint(*_dataHist, _depList, _histObsList, _intOrder, obs, xlo, xhi);
}


////////////////////////////////////////////////////////////////////////////////
/// Return sampling hint for making curves of (projections) of this function
/// as the recursive division strategy of RooCurve cannot deal efficiently
/// with the vertical lines that occur in a non-interpolated histogram
/// Return the histogram bin boundaries mapped into the plot-observable
/// coordinate. Handles the case where the function observable is a
/// transformation (shift, scale, ...) of the plot observable so that the
/// returned boundaries are placed at the correct positions in the plot.

std::list<double>* RooHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
{
Expand All @@ -339,86 +339,21 @@ std::list<double>* RooHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo,
return nullptr ;
}

// Find histogram observable corresponding to pdf observable
RooAbsArg* hobs(nullptr) ;
for (auto i = 0u; i < _histObsList.size(); ++i) {
const auto harg = _histObsList[i];
const auto parg = _depList[i];
if (std::string(parg->GetName())==obs.GetName()) {
hobs=harg ;
}
const std::vector<double> boundaries =
RooHistPdf::histogramBoundariesInPlotObs(*_dataHist, _depList, _histObsList, obs, xlo, xhi);
if (boundaries.empty()) {
return nullptr;
}

// std::cout << "RooHistFunc::bb(" << GetName() << ") histObs = " << _histObsList << std::endl ;
// std::cout << "RooHistFunc::bb(" << GetName() << ") pdfObs = " << _depList << std::endl ;

RooAbsRealLValue* transform = nullptr;
if (!hobs) {

// Considering alternate: input observable is histogram observable and pdf observable is transformation in terms of it
RooAbsArg* pobs = nullptr;
for (auto i = 0u; i < _histObsList.size(); ++i) {
const auto harg = _histObsList[i];
const auto parg = _depList[i];
if (std::string(harg->GetName())==obs.GetName()) {
pobs=parg ;
hobs=harg ;
}
}

// Not found, or check that matching pdf observable is an l-value dependent on histogram observable fails
if (!hobs || !(pobs->dependsOn(obs) && dynamic_cast<RooAbsRealLValue*>(pobs))) {
std::cout << "RooHistFunc::binBoundaries(" << GetName() << ") obs = " << obs.GetName() << " hobs is not found, returning null" << std::endl ;
return nullptr ;
}

// Now we are in business - we are in a situation where the pdf observable LV(x), mapping to a histogram observable x
// We can return bin boundaries by mapping the histogram boundaties through the inverse of the LV(x) transformation
transform = dynamic_cast<RooAbsRealLValue*>(pobs) ;
}


// std::cout << "hobs = " << hobs->GetName() << std::endl ;
// std::cout << "transform = " << (transform?transform->GetName():"<none>") << std::endl ;

// Check that observable is in dataset, if not no hint is generated
RooAbsArg* xtmp = _dataHist->get()->find(hobs->GetName()) ;
if (!xtmp) {
std::cout << "RooHistFunc::binBoundaries(" << GetName() << ") hobs = " << hobs->GetName() << " is not found in dataset?" << std::endl ;
_dataHist->get()->Print("v") ;
return nullptr ;
}
RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dataHist->get()->find(hobs->GetName())) ;
if (!lvarg) {
std::cout << "RooHistFunc::binBoundaries(" << GetName() << ") hobs = " << hobs->GetName() << " but is not an LV, returning null" << std::endl ;
return nullptr ;
}

// Retrieve position of all bin boundaries
const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
double* boundaries = binning->array() ;

auto hint = new std::list<double> ;

double delta = (xhi-xlo)*1e-8 ;

// Construct array with pairs of points positioned epsilon to the left and
// right of the bin boundaries
for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
if (boundaries[i]>xlo-delta && boundaries[i]<xhi+delta) {

double boundary = boundaries[i] ;
if (transform) {
transform->setVal(boundary) ;
//cout << "transform bound " << boundary << " using " << transform->GetName() << " result " << obs.getVal() << std::endl ;
hint->push_back(obs.getVal()) ;
} else {
hint->push_back(boundary) ;
}
}
auto hint = new std::list<double>;
const double delta = (xhi - xlo) * 1e-8;
for (double b : boundaries) {
if (b > xlo - delta && b < xhi + delta) {
hint->push_back(b);
}
}

return hint ;
return hint;
}


Expand Down
151 changes: 108 additions & 43 deletions roofit/roofitcore/src/RooHistPdf.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -416,6 +416,97 @@ bool RooHistPdf::forceAnalyticalInt(const RooAbsArg& dep) const
return forceAnalyticalInt(_pdfObsList, dep);
}

////////////////////////////////////////////////////////////////////////////////
/// Return the histogram bin boundaries mapped to the coordinates of the plot
/// observable `obs`. This handles the case where the p.d.f. observable is a
/// transformation of the plot observable (for example a shift like
/// `x_shifted = x - shift`), so that the returned boundaries are the positions
/// of the histogram bin edges as they appear when the p.d.f. is plotted as a
/// function of `obs`.
///
/// Returns an empty vector if no corresponding histogram observable can be
/// located, or if a non-invertible transformation is detected.

std::vector<double> RooHistPdf::histogramBoundariesInPlotObs(RooDataHist const &dataHist, RooArgSet const &pdfObsList,
RooArgSet const &histObsList, RooAbsRealLValue &obs,
double xlo, double xhi)
{
std::vector<double> result;

// Locate the histogram observable that corresponds to the plot observable.
// These are the scenarios we handle:
// 1. The plot observable `obs` has the same name as one of the p.d.f.
// observables. That p.d.f. observable is the histogram observable
// itself (no transformation).
// 2. The p.d.f. observable is a non-trivial function of `obs` (e.g. a
// shift or scale). We have to invert that function to map the
// histogram bin boundaries back into the plot-observable coordinate.
RooAbsArg *histObs = nullptr;
RooAbsArg *pdfObs = nullptr;
bool useTransform = false;
for (unsigned int i = 0; i < pdfObsList.size(); ++i) {
if (std::string(obs.GetName()) == pdfObsList[i]->GetName()) {
histObs = histObsList[i];
break;
}
}
if (!histObs) {
for (unsigned int i = 0; i < pdfObsList.size(); ++i) {
auto *po = dynamic_cast<RooAbsReal *>(pdfObsList[i]);
if (po && po->dependsOn(obs)) {
histObs = histObsList[i];
pdfObs = pdfObsList[i];
useTransform = true;
break;
}
}
}
if (!histObs) {
return result;
}

auto *lval = dynamic_cast<RooAbsLValue *>(dataHist.get()->find(histObs->GetName()));
if (!lval) {
return result;
}

const RooAbsBinning *binning = lval->getBinningPtr(nullptr);
const int nBoundaries = binning->numBoundaries();
const double *rawBoundaries = binning->array();

if (!useTransform) {
result.assign(rawBoundaries, rawBoundaries + nBoundaries);
return result;
}

// Invert the transformation `pdfObs = T(obs)` by linearly approximating
// it on the range of `obs` and solving for `obs` at each histogram bin
// boundary. This handles both l-values such as RooLinearVar and generic
// RooAbsReals such as RooFormulaVars that describe a shift or scale
// (the approximation is exact for linear transformations).
auto *pdfObsReal = static_cast<RooAbsReal *>(pdfObs);
const double probeLo = std::max(xlo, obs.getMin());
const double probeHi = std::min(xhi, obs.getMax());
if (probeHi > probeLo) {
const double obsSave = obs.getVal();
obs.setVal(probeLo);
const double plo = pdfObsReal->getVal();
obs.setVal(probeHi);
const double phi = pdfObsReal->getVal();
obs.setVal(obsSave);
const double slope = (phi - plo) / (probeHi - probeLo);
if (slope != 0.0) {
result.reserve(nBoundaries);
for (int i = 0; i < nBoundaries; ++i) {
result.push_back(probeLo + (rawBoundaries[i] - plo) / slope);
}
}
}

// The transformation may reverse the order (negative slope).
std::sort(result.begin(), result.end());
return result;
}

////////////////////////////////////////////////////////////////////////////////
/// Return sampling hint for making curves of (projections) of this function
Expand All @@ -441,40 +532,21 @@ std::list<double>* RooHistPdf::plotSamplingHint(RooDataHist const& dataHist,
return nullptr;
}

// Check that observable is in dataset, if not no hint is generated
RooAbsArg* dhObs = nullptr;
for (unsigned int i=0; i < pdfObsList.size(); ++i) {
RooAbsArg* histObs = histObsList[i];
RooAbsArg* pdfObs = pdfObsList[i];
if (std::string(obs.GetName())==pdfObs->GetName()) {
dhObs = dataHist.get()->find(histObs->GetName()) ;
break;
}
}

if (!dhObs) {
return nullptr;
}
RooAbsLValue* lval = dynamic_cast<RooAbsLValue*>(dhObs) ;
if (!lval) {
return nullptr;
const std::vector<double> boundaries = histogramBoundariesInPlotObs(dataHist, pdfObsList, histObsList, obs, xlo, xhi);
if (boundaries.empty()) {
return nullptr;
}

// Retrieve position of all bin boundaries

const RooAbsBinning* binning = lval->getBinningPtr(nullptr);
std::span<const double> boundaries{binning->array(), static_cast<std::size_t>(binning->numBoundaries())};

// 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 sampling hint for making curves of (projections) of this function
/// as the recursive division strategy of RooCurve cannot deal efficiently
/// with the vertical lines that occur in a non-interpolated histogram
/// Return the histogram bin boundaries in the plot-observable coordinate.
/// The returned list is used by RooCurve's adaptive sampler to avoid stepping
/// over the discontinuities of a non-interpolated histogram, and by
/// RooBinIntegrator when integrating over this p.d.f.

std::list<double>* RooHistPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
{
Expand All @@ -483,27 +555,20 @@ std::list<double>* RooHistPdf::binBoundaries(RooAbsRealLValue& obs, double xlo,
return nullptr;
}

// Check that observable is in dataset, if not no hint is generated
RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dataHist->get()->find(obs.GetName())) ;
if (!lvarg) {
return nullptr ;
const std::vector<double> boundaries =
histogramBoundariesInPlotObs(*_dataHist, _pdfObsList, _histObsList, obs, xlo, xhi);
if (boundaries.empty()) {
return nullptr;
}

// Retrieve position of all bin boundaries
const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
double* boundaries = binning->array() ;

auto hint = new std::list<double> ;

// Construct array with pairs of points positioned epsilon to the left and
// right of the bin boundaries
for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
if (boundaries[i]>=xlo && boundaries[i]<=xhi) {
hint->push_back(boundaries[i]) ;
}
auto hint = new std::list<double>;
for (double b : boundaries) {
if (b >= xlo && b <= xhi) {
hint->push_back(b);
}
}

return hint ;
return hint;
}


Expand Down
Loading
Loading