RooBinSamplingPdf.cxx

Go to the documentation of this file.

// Authors: Stephan Hageboeck, CERN; Andrea Sciandra, SCIPP-UCSC/Atlas; Nov 2020
 
/*****************************************************************************
 * RooFit
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2020, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/
 
/**
 * \class RooBinSamplingPdf
 * The RooBinSamplingPdf is supposed to be used as an adapter between a continuous PDF
 * and a binned distribution.
 * When RooFit is used to fit binned data, and the PDF is continuous, it takes the probability density
 * at the bin centre as a proxy for the probability averaged (integrated) over the entire bin. This is
 * correct only if the second derivative of the function vanishes, though. This is shown in the plots
 * below.
 *
 * For PDFs that have larger curvatures, the RooBinSamplingPdf can be used. It integrates the PDF in each
 * bin using an adaptive integrator. This usually requires 21 times more function evaluations, but significantly
 * reduces biases due to better sampling of the PDF. The integrator can be accessed from the outside
 * using integrator(). This can be used to change the integration rules, so less/more function evaluations are
 * performed. The target precision of the integrator can be set in the constructor.
 *
 * If the wrapped PDF supports analytical integration over the sampled observable, the exact analytical integral
 * is used for each bin instead of the numeric integrator. This is both faster and more accurate, and happens
 * transparently without any user intervention.
 *
 *
 * ### How to use it
 * There are two ways to use this class:
 * - Manually wrap a PDF:
 * ```
 *   RooBinSamplingPdf binSampler("<name>", "title", <binned observable of PDF>, <original PDF> [, <precision for integrator>]);
 *   binSampler.fitTo(data);
 * ```
 *   When a PDF is wrapped with a RooBinSamplingPDF, just use the bin sampling PDF instead of the original one for fits
 *   or plotting etc.
 *   \note The binning will be taken from the observable. Make sure that this binning is the same as the one of the dataset that should be fit.
 *   Use RooRealVar::setBinning() to adapt it.
 * - Instruct test statistics to carry out this wrapping automatically:
 * ```
 *   pdf.fitTo(data, IntegrateBins(<precision>));
 * ```
 *   This method is especially useful when used with a simultaneous PDF, since each component will automatically be wrapped,
 *   depending on the value of `precision`:
 *   - `precision < 0.`: None of the PDFs are touched, bin sampling is off.
 *   - `precision = 0.`: Continuous PDFs that are fit to a RooDataHist are wrapped into a RooBinSamplingPdf. The target precision
 *      forwarded to the integrator is 1.E-4 (the default argument of the constructor).
 *   - `precision > 0.`: All continuous PDFs are automatically wrapped into a RooBinSamplingPdf, regardless of what data they are
 *      fit to (see next paragraph). The same `'precision'` is used for all integrators.
 *
 * ### Simulating a binned fit using RooDataSet
 *   Some frameworks use unbinned data (RooDataSet) to simulate binned datasets. By adding one entry for each bin centre with the
 *   appropriate weight, one can achieve the same result as fitting with RooDataHist. In this case, however, RooFit cannot
 *   auto-detect that a binned fit is running, and that an integration over the bin is desired (note that there are no bins to
 *   integrate over in this kind of dataset).
 *
 *   In this case, `IntegrateBins(>0.)` needs to be used, and the desired binning needs to be assigned to the observable
 *   of the dataset:
 *   ```
 *     RooRealVar x("x", "x", 0., 5.);
 *     x.setBins(10);
 *
 *     // <create dataset and model>
 *
 *     model.fitTo(data, IntegrateBins(>0.));
 *   ```
 *
 *   \see RooAbsPdf::fitTo()
 *   \see IntegrateBins()
 *
 *   \note This feature is currently limited to one-dimensional PDFs.
 *
 *
 * \htmlonly <style>div.image img[src="RooBinSamplingPdf_OFF.png"]{width:12cm;}</style> \endhtmlonly
 * \htmlonly <style>div.image img[src="RooBinSamplingPdf_ON.png" ]{width:12cm;}</style> \endhtmlonly
 * <table>
 * <tr><th> Binned fit without %RooBinSamplingPdf    <th> Binned fit with %RooBinSamplingPdf   </td></tr>
 * <tr><td> \image html RooBinSamplingPdf_OFF.png ""
 *   </td>
 * <td> \image html RooBinSamplingPdf_ON.png ""
 *   </td></tr>
 * </table>
 *
 */
 
 
#include "RooBinSamplingPdf.h"
 
#include "RooFitImplHelpers.h"
#include "RooRealBinding.h"
#include "RooRealVar.h"
#include "RooGlobalFunc.h"
#include "RooDataHist.h"
 
#include "Math/Integrator.h"
 
#include <algorithm>
 
////////////////////////////////////////////////////////////////////////////////
/// Construct a new RooBinSamplingPdf.
/// \param[in] name A name to identify this object.
/// \param[in] title Title (for e.g. plotting)
/// \param[in] observable Observable to integrate over (the one that is binned).
/// \param[in] inputPdf A PDF whose bins should be sampled with higher precision.
/// \param[in] epsilon Relative precision for the integrator, which is used to sample the bins.
/// Note that ROOT's default is to use an adaptive integrator, which in its first iteration usually reaches
/// relative precision of 1.E-4 or better. Therefore, asking for lower precision rarely has an effect.
RooBinSamplingPdf::RooBinSamplingPdf(const char *name, const char *title, RooAbsRealLValue& observable,
    RooAbsPdf& inputPdf, double epsilon) :
      RooAbsPdf(name, title),
      _pdf("inputPdf", "Function to be converted into a PDF", this, inputPdf),
      _observable("observable", "Observable to integrate over", this, observable, true, true),
      _relEpsilon(epsilon) {
  if (!_pdf->dependsOn(*_observable)) {
    throw std::invalid_argument(std::string("RooBinSamplingPDF(") + GetName()
        + "): The PDF " + _pdf->GetName() + " needs to depend on the observable "
        + _observable->GetName());
  }
}
 
 
 ////////////////////////////////////////////////////////////////////////////////
 /// Copy a RooBinSamplingPdf.
 /// \param[in] other PDF to copy.
 /// \param[in] name Optionally rename the copy.
 RooBinSamplingPdf::RooBinSamplingPdf(const RooBinSamplingPdf& other, const char* name) :
   RooAbsPdf(other, name),
   _pdf("inputPdf", this, other._pdf),
   _observable("observable", this, other._observable),
   _relEpsilon(other._relEpsilon) { }
 
 
////////////////////////////////////////////////////////////////////////////////
/// Integrate the PDF over the current bin of the observable.
double RooBinSamplingPdf::evaluate() const {
  const unsigned int bin = _observable->getBin();
  const double low = _observable->getBinning().binLow(bin);
  const double high = _observable->getBinning().binHigh(bin);
 
  const double oldX = _observable->getVal();
  double result;
  {
    // Important: When the integrator samples x, caching of sub-tree values needs to be off.
    DisableCachingRAII disableCaching(inhibitDirty());
    result = integrate(_normSet, low, high) / (high-low);
  }
 
  _observable->setVal(oldX);
 
  return result;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Integrate the PDF over all its bins, and return a batch with those values.
/// \param[in,out] ctx Struct with evaluation data
void RooBinSamplingPdf::doEval(RooFit::EvalContext &ctx) const
{
   std::span<double> output = ctx.output();
 
   // Retrieve binning, which we need to compute the probabilities
   auto boundaries = binBoundaries();
   auto xValues = ctx.at(_observable);
 
   // Important: When the integrator samples x, caching of sub-tree values needs to be off.
   DisableCachingRAII disableCaching(inhibitDirty());
 
   // Now integrate PDF in each bin:
   for (unsigned int i = 0; i < xValues.size(); ++i) {
      const double x = xValues[i];
      const auto upperIt = std::upper_bound(boundaries.begin(), boundaries.end(), x);
      const unsigned int bin = std::distance(boundaries.begin(), upperIt) - 1;
      assert(bin < boundaries.size());
 
      output[i] = integrate(nullptr, boundaries[bin], boundaries[bin + 1]) / (boundaries[bin + 1] - boundaries[bin]);
   }
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Get the bin boundaries for the observable.
/// These will be recomputed whenever the shape of this object is dirty.
std::span<const double> RooBinSamplingPdf::binBoundaries() const {
  if (isShapeDirty() || _binBoundaries.empty()) {
    _binBoundaries.clear();
    const RooAbsBinning& binning = _observable->getBinning(nullptr);
    const double* boundaries = binning.array();
 
    for (int i=0; i < binning.numBoundaries(); ++i) {
      _binBoundaries.push_back(boundaries[i]);
    }
 
    assert(std::is_sorted(_binBoundaries.begin(), _binBoundaries.end()));
 
    clearShapeDirty();
  }
 
  return {_binBoundaries};
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Return a list of all bin boundaries, so the PDF is plotted correctly.
/// \param[in] obs Observable to generate the boundaries for.
/// \param[in] xlo Beginning of range to create list of boundaries for.
/// \param[in] xhi End of range to create to create list of boundaries for.
/// \return Pointer to a list to be deleted by caller.
std::list<double>* RooBinSamplingPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const {
  if (obs.namePtr() != _observable->namePtr()) {
    coutE(Plotting) << "RooBinSamplingPdf::binBoundaries(" << GetName() << "): observable '" << obs.GetName()
        << "' is not the observable of this PDF ('" << _observable->GetName() << "')." << std::endl;
    return nullptr;
  }
 
  auto list = new std::list<double>;
  for (double val : binBoundaries()) {
    if (xlo <= val && val < xhi)
      list->push_back(val);
  }
 
  return list;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Return a list of all bin edges, so the PDF is plotted as a step function.
/// \param[in] obs Observable to generate the sampling hint for.
/// \param[in] xlo Beginning of range to create sampling hint for.
/// \param[in] xhi End of range to create sampling hint for.
/// \return Pointer to a list to be deleted by caller.
std::list<double>* RooBinSamplingPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const {
  if (obs.namePtr() != _observable->namePtr()) {
    coutE(Plotting) << "RooBinSamplingPdf::plotSamplingHint(" << GetName() << "): observable '" << obs.GetName()
        << "' is not the observable of this PDF ('" << _observable->GetName() << "')." << std::endl;
    return nullptr;
  }
 
  auto binEdges = new std::list<double>;
  const auto& binning = obs.getBinning();
 
  for (unsigned int bin=0, numBins = static_cast<unsigned int>(binning.numBins()); bin < numBins; ++bin) {
    const double low  = std::max(binning.binLow(bin), xlo);
    const double high = std::min(binning.binHigh(bin), xhi);
    const double width = high - low;
 
    // Check if this bin is in plotting range at all
    if (low >= high)
      continue;
 
    // Move support points slightly inside the bin, so step function is plotted correctly.
    binEdges->push_back(low  + 0.001 * width);
    binEdges->push_back(high - 0.001 * width);
  }
 
  return binEdges;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Direct access to the unique_ptr holding the integrator that's used to sample the bins.
/// This can be used to change options such as sampling accuracy or to entirely exchange the integrator.
///
/// #### Example: Use the 61-point Gauss-Kronrod integration rule
/// ```{.cpp}
///   ROOT::Math::IntegratorOneDimOptions intOptions = pdf.integrator()->Options();
///   intOptions.SetNPoints(6); // 61-point integration rule
///   intOptions.SetRelTolerance(1.E-9); // Smaller tolerance -> more subdivisions
///   pdf.integrator()->SetOptions(intOptions);
/// ```
/// \note see ROOT::Math::IntegratorOneDim::SetOptions for more details on integration options.
/// \note When RooBinSamplingPdf is loaded from files, integrator options will fall back to the default values.
std::unique_ptr<ROOT::Math::IntegratorOneDim>& RooBinSamplingPdf::integrator() const {
  if (!_integrator) {
    _integrator = std::make_unique<ROOT::Math::IntegratorOneDim>(*this,
        ROOT::Math::IntegrationOneDim::kADAPTIVE, // GSL Integrator. Will really get it only if MathMore enabled.
        -1., _relEpsilon, // Abs epsilon = default, rel epsilon set by us.
        0, // We don't limit the sub-intervals. Steer run time via _relEpsilon.
        2 // This should read ROOT::Math::Integration::kGAUSS21, but this is in MathMore, so we cannot include it here.
        );
  }
 
  return _integrator;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Binding used by the integrator to evaluate the PDF.
double RooBinSamplingPdf::operator()(double x) const {
  _observable->setVal(x);
  return _pdf;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Check once whether the wrapped PDF can integrate over the observable
/// analytically. If so, the analytical integral code is cached so that
/// integrate() can use the exact integral instead of the numeric integrator.
void RooBinSamplingPdf::initializeAnalyticalIntegral(const RooArgSet* normSet) const {
  // Setting a named range to select the bin boundaries requires a RooRealVar.
  // For other observable types we stick to the numeric integrator.
  auto *observable = dynamic_cast<RooRealVar *>(&*_observable);
  if (!observable) {
    _analyticalIntegralCode = 0;
    return;
  }
 
  _analyticalIntegralRangeName = std::string("_binSampling_") + GetName();
 
  // Define the range so that PDFs that inspect it in getAnalyticalIntegral()
  // find a valid one. The actual bin boundaries are filled in for each bin in
  // integrate().
  observable->setRange(_analyticalIntegralRangeName.c_str(), observable->getMin(), observable->getMax());
 
  RooArgSet allVars{*_observable};
  RooArgSet analVars;
  _analyticalIntegralCode = _pdf->getAnalyticalIntegralWN(allVars, analVars, normSet, _analyticalIntegralRangeName.c_str());
 
  // Only use analytical integration if the observable is really integrated
  // analytically. Otherwise, fall back to the numeric integrator.
  if (_analyticalIntegralCode != 0 && !analVars.contains(*_observable)) {
    _analyticalIntegralCode = 0;
  }
 
  if (_analyticalIntegralCode != 0) {
    coutI(NumIntegration) << "RooBinSamplingPdf::integrate(" << GetName()
        << "): using the analytical integral of " << _pdf->GetName()
        << " to sample the bins instead of the numeric integrator." << std::endl;
  }
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Integrate the wrapped PDF over a single bin, with the given norm set and limits.
/// If the wrapped PDF supports analytical integration over the observable, the
/// exact integral is used. Otherwise, the numeric integrator is employed.
///
/// The result must match the numeric path, which integrates the value returned by
/// operator(), i.e. the PDF value normalized over `normSet`. Therefore, the
/// analytical integral is normalized over `normSet` as well via analyticalIntegralWN().
double RooBinSamplingPdf::integrate(const RooArgSet* normSet, double low, double high) const {
  if (_analyticalIntegralCode == -1) {
    initializeAnalyticalIntegral(normSet);
  }
 
  if (_analyticalIntegralCode != 0) {
    // The analytical path is only enabled for RooRealVar observables (see
    // initializeAnalyticalIntegral()), so this static_cast is safe.
    static_cast<RooRealVar &>(*_observable).setRange(_analyticalIntegralRangeName.c_str(), low, high);
    return _pdf->analyticalIntegralWN(_analyticalIntegralCode, normSet, _analyticalIntegralRangeName.c_str());
  }
 
  return integrator()->Integral(low, high);
}
 
 
/// Creates a wrapping RooBinSamplingPdf if appropriate.
/// \param[in] pdf The input pdf.
/// \param[in] data The dataset to be used in the fit, used to figure out the
///            observables and whether the dataset is binned.
/// \param[in] precision Precision argument for all created RooBinSamplingPdfs.
std::unique_ptr<RooAbsPdf> RooBinSamplingPdf::create(RooAbsPdf& pdf, RooAbsData const &data, double precision) {
  if (precision < 0.)
    return nullptr;
 
  std::unique_ptr<RooArgSet> funcObservables( pdf.getObservables(data) );
  const bool oneDimAndBinned = (1 == std::count_if(funcObservables->begin(), funcObservables->end(), [](const RooAbsArg* arg) {
    auto var = dynamic_cast<const RooRealVar*>(arg);
    return var && var->numBins() > 1;
  }));
 
  if (!oneDimAndBinned) {
    if (precision > 0.) {
      oocoutE(&pdf, Fitting)
          << "Integration over bins was requested, but this is currently only implemented for 1-D fits." << std::endl;
    }
    return nullptr;
  }
 
  // Find the real-valued observable. We don't care about categories.
  auto theObs = std::find_if(funcObservables->begin(), funcObservables->end(), [](const RooAbsArg* arg){
    return dynamic_cast<const RooAbsRealLValue*>(arg);
  });
  assert(theObs != funcObservables->end());
 
  std::unique_ptr<RooAbsPdf> newPdf;
 
  if (precision > 0.) {
    // User forced integration. Let just apply it.
    newPdf = std::make_unique<RooBinSamplingPdf>(
        (std::string(pdf.GetName()) + "_binSampling").c_str(), pdf.GetTitle(),
        *static_cast<RooAbsRealLValue *>(*theObs), pdf, precision);
  } else if (dynamic_cast<RooDataHist const *>(&data) != nullptr &&
             precision == 0. && !pdf.isBinnedDistribution(*data.get())) {
    // User didn't forbid integration, and it seems appropriate with a
    // RooDataHist.
    newPdf = std::make_unique<RooBinSamplingPdf>(
        (std::string(pdf.GetName()) + "_binSampling").c_str(), pdf.GetTitle(),
        *static_cast<RooAbsRealLValue *>(*theObs), pdf);
  }
 
  return newPdf;
}