RooAbsOptTestStatistic.cxx

Go to the documentation of this file.

/// \cond ROOFIT_INTERNAL
 
/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, 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)             *
 *****************************************************************************/
 
/**
\file RooAbsOptTestStatistic.cxx
\class RooAbsOptTestStatistic
\ingroup Roofitcore
 
Abstract base class for test
statistics objects that evaluate a function or PDF at each point of a given
dataset.  This class provides generic optimizations, such as
caching and precalculation of constant terms that can be made for
all such quantities.
 
Implementations should define evaluatePartition(), which calculates the
value of a (sub)range of the dataset and optionally combinedValue(),
which combines the values calculated for each partition. If combinedValue()
is not overloaded, the default implementation will add the partition results
to obtain the combined result.
 
Support for calculation in partitions is needed to allow multi-core
parallelized calculation of test statistics.
**/
 
#include "RooAbsOptTestStatistic.h"
 
#include "Riostream.h"
#include "TClass.h"
#include <cstring>
 
#include "RooAbsData.h"
#include "RooAbsDataStore.h"
#include "RooAbsPdf.h"
#include "RooAddPdf.h"
#include "RooArgSet.h"
#include "RooBinSamplingPdf.h"
#include "RooBinning.h"
#include "RooCategory.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "RooErrorHandler.h"
#include "RooFitImplHelpers.h"
#include "RooGlobalFunc.h"
#include "RooMsgService.h"
#include "RooProdPdf.h"
#include "RooProduct.h"
#include "RooRealSumPdf.h"
#include "RooRealVar.h"
#include "RooVectorDataStore.h"
 
#include "ROOT/StringUtils.hxx"
 
using std::ostream;
 
////////////////////////////////////////////////////////////////////////////////
/// Create a test statistic, and optimise its calculation.
/// \param[in] name Name of the instance.
/// \param[in] title Title (for e.g. plotting).
/// \param[in] real Function to evaluate.
/// \param[in] indata Dataset for which to compute test statistic.
/// \param[in] projDeps A set of projected observables.
/// \param[in] cfg the statistic configuration
///
/// cfg contains:
/// - rangeName If not null, only events in the dataset inside the range will be used in the test
///   statistic calculation.
/// - addCoefRangeName If not null, all RooAddPdf components of `real` will be
///   instructed to fix their fraction definitions to the given named range.
/// - nCPU If > 1, the test statistic calculation will be parallelised over multiple processes. By default, the data
///   is split with 'bulk' partitioning (each process calculates a contiguous block of fraction 1/nCPU
///   of the data). For binned data, this approach may be suboptimal as the number of bins with >0 entries
///   in each processing block may vary greatly; thereby distributing the workload rather unevenly.
/// - interleave Strategy how to distribute events among workers. If an interleave partitioning strategy is used where each partition
///   i takes all bins for which (ibin % ncpu == i), an even distribution of work is more likely.
/// - splitCutRange If true, a different rangeName constructed as `rangeName_{catName}` will be used
///   as range definition for each index state of a RooSimultaneous.
/// - cloneInputData Not used. Data is always cloned.
/// - integrateOverBinsPrecision If > 0, PDF in binned fits are integrated over the bins. This sets the precision. If = 0,
///   only unbinned PDFs fit to RooDataHist are integrated. If < 0, PDFs are never integrated.
RooAbsOptTestStatistic::RooAbsOptTestStatistic(const char *name, const char *title, RooAbsReal &real,
                                               RooAbsData &indata, const RooArgSet &projDeps,
                                               RooAbsTestStatistic::Configuration const &cfg)
   : RooAbsTestStatistic(name, title, real, indata, projDeps, cfg),
     _integrateBinsPrecision(cfg.integrateOverBinsPrecision)
{
   // Don't do a thing in master mode
   if (operMode() != Slave) {
      return;
   }
 
   initSlave(real, indata, projDeps, _rangeName.c_str(), _addCoefRangeName.c_str());
}
 
////////////////////////////////////////////////////////////////////////////////
/// Copy constructor
 
RooAbsOptTestStatistic::RooAbsOptTestStatistic(const RooAbsOptTestStatistic &other, const char *name)
   : RooAbsTestStatistic(other, name),
     _sealed(other._sealed),
     _sealNotice(other._sealNotice),
     _integrateBinsPrecision(other._integrateBinsPrecision)
{
   // Don't do a thing in master mode
   if (operMode() != Slave) {
 
      if (other._normSet) {
         _normSet = new RooArgSet;
         other._normSet->snapshot(*_normSet);
      }
      return;
   }
 
   initSlave(*other._funcClone, *other._dataClone, other._projDeps ? *other._projDeps : RooArgSet(),
             other._rangeName.c_str(), other._addCoefRangeName.c_str());
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
 
void RooAbsOptTestStatistic::initSlave(RooAbsReal& real, RooAbsData& indata, const RooArgSet& projDeps, const char* rangeName,
                   const char* addCoefRangeName) {
  // ******************************************************************
  // *** PART 1 *** Clone incoming pdf, attach to each other *
  // ******************************************************************
 
  // Clone FUNC
  _funcClone = RooHelpers::cloneTreeWithSameParameters(real, indata.get()).release();
  _funcCloneSet = nullptr ;
 
  // Attach FUNC to data set
  _funcObsSet = std::unique_ptr<RooArgSet>{_funcClone->getObservables(indata)}.release();
 
  if (_funcClone->getAttribute("BinnedLikelihood")) {
    _funcClone->setAttribute("BinnedLikelihoodActive") ;
  }
 
  // Mark all projected dependents as such
  if (!projDeps.empty()) {
    std::unique_ptr<RooArgSet> projDataDeps{_funcObsSet->selectCommon(projDeps)};
    projDataDeps->setAttribAll("projectedDependent") ;
  }
 
  // If PDF is a RooProdPdf (with possible constraint terms)
  // analyze pdf for actual parameters (i.e those in unconnected constraint terms should be
  // ignored as here so that the test statistic will not be recalculated if those
  // are changed
  RooProdPdf* pdfWithCons = dynamic_cast<RooProdPdf*>(_funcClone) ;
  if (pdfWithCons) {
 
    std::unique_ptr<RooArgSet> connPars{pdfWithCons->getConnectedParameters(*indata.get())};
    // Add connected parameters as servers
    _paramSet.add(*connPars) ;
 
  } else {
    // Add parameters as servers
    _funcClone->getParameters(indata.get(), _paramSet);
  }
 
  // Store normalization set
  _normSet = new RooArgSet;
  indata.get()->snapshot(*_normSet, false);
 
  // Expand list of observables with any observables used in parameterized ranges.
  // This NEEDS to be a counting loop since we are inserting during the loop.
  for (std::size_t i = 0; i < _funcObsSet->size(); ++i) {
    auto realDepRLV = dynamic_cast<const RooAbsRealLValue*>((*_funcObsSet)[i]);
    if (realDepRLV && realDepRLV->isDerived()) {
      RooArgSet tmp2;
      realDepRLV->leafNodeServerList(&tmp2, nullptr, true);
      _funcObsSet->add(tmp2,true);
    }
  }
 
 
 
  // ******************************************************************
  // *** PART 2 *** Clone and adjust incoming data, attach to PDF     *
  // ******************************************************************
 
  // Check if the fit ranges of the dependents in the data and in the FUNC are consistent
  const RooArgSet* dataDepSet = indata.get() ;
  for (const auto arg : *_funcObsSet) {
 
    // Check that both dataset and function argument are of type RooRealVar
    RooRealVar* realReal = dynamic_cast<RooRealVar*>(arg) ;
    if (!realReal) continue ;
    RooRealVar* datReal = dynamic_cast<RooRealVar*>(dataDepSet->find(realReal->GetName())) ;
    if (!datReal) continue ;
 
    // Check that range of observables in pdf is equal or contained in range of observables in data
 
    if (!realReal->getBinning().lowBoundFunc() && realReal->getMin()<(datReal->getMin()-1e-6)) {
      coutE(InputArguments) << "RooAbsOptTestStatistic: ERROR minimum of FUNC observable " << arg->GetName()
                 << "(" << realReal->getMin() << ") is smaller than that of "
                 << arg->GetName() << " in the dataset (" << datReal->getMin() << ")" << std::endl ;
      RooErrorHandler::softAbort() ;
      return ;
    }
 
    if (!realReal->getBinning().highBoundFunc() && realReal->getMax()>(datReal->getMax()+1e-6)) {
      coutE(InputArguments) << "RooAbsOptTestStatistic: ERROR maximum of FUNC observable " << arg->GetName()
                 << " is larger than that of " << arg->GetName() << " in the dataset" << std::endl ;
      RooErrorHandler::softAbort() ;
      return ;
    }
  }
 
  // Copy data and strip entries lost by adjusted fit range, _dataClone ranges will be copied from realDepSet ranges
  if (rangeName && strlen(rangeName)) {
    _dataClone = std::unique_ptr<RooAbsData>{indata.reduce(RooFit::SelectVars(*_funcObsSet),RooFit::CutRange(rangeName))}.release();
  } else {
    _dataClone = static_cast<RooAbsData*>(indata.Clone()) ;
  }
  _ownData = true ;
 
 
  // ******************************************************************
  // *** PART 3 *** Make adjustments for fit ranges, if specified     *
  // ******************************************************************
 
  std::unique_ptr<RooArgSet> origObsSet( real.getObservables(indata) );
  if (rangeName && strlen(rangeName)) {
    cxcoutI(Fitting) << "RooAbsOptTestStatistic::ctor(" << GetName() << ") constructing test statistic for sub-range named " << rangeName << std::endl ;
 
    if(auto pdfClone = dynamic_cast<RooAbsPdf*>(_funcClone)) {
       pdfClone->setNormRange(rangeName);
    }
 
    // Print warnings if the requested ranges are not available for the observable
    for (const auto arg : *_funcObsSet) {
 
      if (auto realObs = dynamic_cast<RooRealVar*>(arg)) {
 
        auto tokens = ROOT::Split(rangeName, ",");
        for(std::string const& token : tokens) {
          if(!realObs->hasRange(token.c_str())) {
             std::stringstream errMsg;
             errMsg << "The observable \"" << realObs->GetName() << "\" doesn't define the requested range \""
                    << token << "\". Replacing it with the default range." << std::endl;
             coutI(Fitting) << errMsg.str() << std::endl;
          }
        }
      }
    }
  }
 
 
  // ******************************************************************
  // *** PART 3.2 *** Binned fits                                     *
  // ******************************************************************
 
  setUpBinSampling();
 
 
  // Fix RooAddPdf coefficients to original normalization range
  if (rangeName && strlen(rangeName)) {
 
    // WVE Remove projected dependents from normalization
    _funcClone->fixAddCoefNormalization(*_dataClone->get(),false) ;
 
    if (addCoefRangeName && strlen(addCoefRangeName)) {
      cxcoutI(Fitting) << "RooAbsOptTestStatistic::ctor(" << GetName()
                 << ") fixing interpretation of coefficients of any RooAddPdf component to range " << addCoefRangeName << std::endl ;
      _funcClone->fixAddCoefRange(addCoefRangeName,false) ;
    }
  }
 
 
  // This is deferred from part 2 - but must happen after part 3 - otherwise invalid bins cannot be properly marked in cacheValidEntries
  _dataClone->attachBuffers(*_funcObsSet) ;
  setEventCount(_dataClone->numEntries()) ;
 
 
 
 
  // *********************************************************************
  // *** PART 4 *** Adjust normalization range for projected observables *
  // *********************************************************************
 
  // Remove projected dependents from normalization set
  if (!projDeps.empty()) {
 
    _projDeps = new RooArgSet;
    projDeps.snapshot(*_projDeps, false) ;
 
    //RooArgSet* tobedel = (RooArgSet*) _normSet->selectCommon(*_projDeps) ;
    _normSet->remove(*_projDeps,true,true) ;
 
    // Mark all projected dependents as such
    RooArgSet projDataDeps;
    _funcObsSet->selectCommon(*_projDeps, projDataDeps);
    projDataDeps.setAttribAll("projectedDependent") ;
  }
 
 
  coutI(Optimization) << "RooAbsOptTestStatistic::ctor(" << GetName() << ") optimizing internal clone of p.d.f for likelihood evaluation."
      << "Lazy evaluation and associated change tracking will disabled for all nodes that depend on observables" << std::endl ;
 
 
  // *********************************************************************
  // *** PART 4 *** Finalization and activation of optimization          *
  // *********************************************************************
 
  // Redirect pointers of base class to clone
  _func = _funcClone ;
  _data = _dataClone ;
 
  _funcClone->getVal(_normSet) ;
 
  optimizeCaching() ;
 
  // It would be unusual if the global observables are used in the likelihood
  // outside of the constraint terms, but if they are we have to be consistent
  // and also redirect them to the snapshots in the dataset if appropriate.
  if(_takeGlobalObservablesFromData && _data->getGlobalObservables()) {
    recursiveRedirectServers(*_data->getGlobalObservables()) ;
  }
 
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Destructor
 
RooAbsOptTestStatistic::~RooAbsOptTestStatistic()
{
  if (operMode()==Slave) {
    delete _funcClone ;
    delete _funcObsSet ;
    if (_projDeps) {
      delete _projDeps ;
    }
    if (_ownData) {
      delete _dataClone ;
    }
  }
  delete _normSet ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Method to combined test statistic results calculated into partitions into
/// the global result. This default implementation adds the partition return
/// values
 
double RooAbsOptTestStatistic::combinedValue(RooAbsReal** array, Int_t n) const
{
  // Default implementation returns sum of components
  double sum(0);
  double carry(0);
  for (Int_t i = 0; i < n; ++i) {
    double y = array[i]->getValV();
    carry += reinterpret_cast<RooAbsOptTestStatistic*>(array[i])->getCarry();
    y -= carry;
    const double t = sum + y;
    carry = (t - sum) - y;
    sum = t;
  }
  _evalCarry = carry;
  return sum ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Catch server redirect calls and forward to internal clone of function
 
bool RooAbsOptTestStatistic::redirectServersHook(const RooAbsCollection& newServerList, bool mustReplaceAll, bool nameChange, bool isRecursive)
{
  RooAbsTestStatistic::redirectServersHook(newServerList,mustReplaceAll,nameChange,isRecursive) ;
  if (operMode()!=Slave) return false ;
  bool ret = _funcClone->recursiveRedirectServers(newServerList,false,nameChange) ;
  return ret || RooAbsReal::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursive);
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Catch print hook function and forward to function clone
 
void RooAbsOptTestStatistic::printCompactTreeHook(ostream& os, const char* indent)
{
  RooAbsTestStatistic::printCompactTreeHook(os,indent) ;
  if (operMode()!=Slave) return ;
  TString indent2(indent) ;
  indent2 += "opt >>" ;
  _funcClone->printCompactTree(os,indent2.Data()) ;
  os << indent2 << " dataset clone = " << _dataClone << " first obs = " << _dataClone->get()->first() << std::endl ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// This method changes the value caching logic for all nodes that depends on any of the observables
/// as defined by the given dataset. When evaluating a test statistic constructed from the RooAbsReal
/// with a dataset the observables are guaranteed to change with every call, thus there is no point
/// in tracking these changes which result in a net overhead. Thus for observable-dependent nodes,
/// the evaluation mechanism is changed from being dependent on a 'valueDirty' flag to guaranteed evaluation.
/// On the dataset side, the observables objects are modified to no longer send valueDirty messages
/// to their client
 
void RooAbsOptTestStatistic::optimizeCaching()
{
  // Trigger create of all object caches now in nodes that have deferred object creation
  // so that cache contents can be processed immediately
  _funcClone->getVal(_normSet) ;
 
  // Set value caching mode for all nodes that depend on any of the observables to ADirty
  _funcClone->optimizeCacheMode(*_funcObsSet) ;
 
  // Disable propagation of dirty state flags for observables
  _dataClone->setDirtyProp(false) ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Change dataset that is used to given one. If cloneData is true, a clone of
/// in the input dataset is made.  If the test statistic was constructed with
/// a range specification on the data, the cloneData argument is ignored and
/// the data is always cloned.
bool RooAbsOptTestStatistic::setDataSlave(RooAbsData& indata, bool cloneData, bool ownNewData)
{
 
  if (operMode()==SimMaster) {
    return false ;
  }
 
 
  // If the current dataset is owned, transfer the ownership to unique pointer
  // that will get out of scope at the end of this function. We can't delete it
  // right now, because there might be global observables in the model that
  // first need to be redirected to the new dataset with a later call to
  // RooAbsArg::recursiveRedirectServers.
  std::unique_ptr<RooAbsData> oldOwnedData;
  if (_ownData) {
    oldOwnedData.reset(_dataClone);
    _dataClone = nullptr ;
  }
 
  if (!cloneData && !_rangeName.empty()) {
    coutW(InputArguments) << "RooAbsOptTestStatistic::setData(" << GetName() << ") WARNING: test statistic was constructed with range selection on data, "
          << "ignoring request to _not_ clone the input dataset" << std::endl ;
    cloneData = true ;
  }
 
  if (cloneData) {
    // Cloning input dataset
    _dataClone = std::unique_ptr<RooAbsData>{indata.reduce(RooFit::SelectVars(*indata.get()),RooFit::CutRange(_rangeName.c_str()))}.release();
    _ownData = true ;
 
  } else {
 
    // Taking input dataset
    _dataClone = &indata ;
    _ownData = ownNewData ;
 
  }
 
  // Attach function clone to dataset
  _dataClone->attachBuffers(*_funcObsSet) ;
  _dataClone->setDirtyProp(false) ;
  _data = _dataClone ;
 
  // Adjust internal event count
  setEventCount(indata.numEntries()) ;
 
  setValueDirty() ;
 
  // It would be unusual if the global observables are used in the likelihood
  // outside of the constraint terms, but if they are we have to be consistent
  // and also redirect them to the snapshots in the dataset if appropriate.
  if(_takeGlobalObservablesFromData && _data->getGlobalObservables()) {
    recursiveRedirectServers(*_data->getGlobalObservables()) ;
  }
 
  return true ;
}
 
 
 
 
////////////////////////////////////////////////////////////////////////////////
 
RooAbsData& RooAbsOptTestStatistic::data()
{
  if (_sealed) {
    bool notice = (sealNotice() && strlen(sealNotice())) ;
    coutW(ObjectHandling) << "RooAbsOptTestStatistic::data(" << GetName()
           << ") WARNING: object sealed by creator - access to data is not permitted: "
           << (notice?sealNotice():"<no user notice>") << std::endl ;
    static RooDataSet dummy ("dummy","dummy",RooArgSet()) ;
    return dummy ;
  }
  return *_dataClone ;
}
 
 
////////////////////////////////////////////////////////////////////////////////
 
const RooAbsData& RooAbsOptTestStatistic::data() const
{
  if (_sealed) {
    bool notice = (sealNotice() && strlen(sealNotice())) ;
    coutW(ObjectHandling) << "RooAbsOptTestStatistic::data(" << GetName()
           << ") WARNING: object sealed by creator - access to data is not permitted: "
           << (notice?sealNotice():"<no user notice>") << std::endl ;
    static RooDataSet dummy ("dummy","dummy",RooArgSet()) ;
    return dummy ;
  }
  return *_dataClone ;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Inspect PDF to find out if we are doing a binned fit to a 1-dimensional unbinned PDF.
/// If this is the case, enable finer sampling of bins by wrapping PDF into a RooBinSamplingPdf.
/// The member _integrateBinsPrecision decides how we act:
/// - < 0: Don't do anything.
/// - = 0: Only enable feature if fitting unbinned PDF to RooDataHist.
/// - > 0: Enable as requested.
void RooAbsOptTestStatistic::setUpBinSampling() {
 
  auto& pdf = static_cast<RooAbsPdf&>(*_funcClone);
  if (auto newPdf = RooBinSamplingPdf::create(pdf, *_dataClone, _integrateBinsPrecision)) {
    newPdf->addOwnedComponents(*_funcClone);
    _funcClone = newPdf.release();
  }
 
}
 
 
/// Returns a suffix string that is unique for RooAbsOptTestStatistic
/// instances that don't share the same cloned input data object.
const char* RooAbsOptTestStatistic::cacheUniqueSuffix() const {
   return Form("_%lx", _dataClone->uniqueId().value()) ;
}
 
/// \endcond