RooAbsL.cxx

Go to the documentation of this file.

/*
 * Project: RooFit
 * Authors:
 *   PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl
 *
 * Copyright (c) 2021, CERN
 *
 * 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)
 */
 
#include <RooFit/TestStatistics/RooAbsL.h>
#include "ConstantTermsOptimizer.h"
#include "RooAbsData.h"
#include "RooAbsDataStore.h" // complete class for access in constOptimizeTestStatistic
 
// for dynamic casts in init_clones:
#include "RooAbsRealLValue.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
 
// other stuff in init_clones:
#include "RooErrorHandler.h"
#include "RooMsgService.h"
 
// concrete classes in getParameters (testing, remove later)
#include "RooRealSumPdf.h"
#include "RooProdPdf.h"
 
namespace RooFit {
namespace TestStatistics {
 
// static function
bool RooAbsL::isExtendedHelper(RooAbsPdf *pdf, Extended extended)
{
   switch (extended) {
   case RooAbsL::Extended::No: {
      return false;
   }
   case RooAbsL::Extended::Yes: {
      return true;
   }
   case RooAbsL::Extended::Auto: {
      return ((pdf->extendMode() == RooAbsPdf::CanBeExtended || pdf->extendMode() == RooAbsPdf::MustBeExtended));
   }
   default: {
      throw std::logic_error("RooAbsL::isExtendedHelper got an unknown extended value!");
   }
   }
}
 
/// After handling cloning (or not) of the pdf and dataset, the public constructors call this private constructor to
/// handle common tasks.
RooAbsL::RooAbsL(std::shared_ptr<RooAbsPdf> pdf, std::shared_ptr<RooAbsData> data, std::size_t N_events,
                 std::size_t N_components, Extended extended)
   : pdf_(std::move(pdf)),
     data_(std::move(data)),
     N_events_(N_events),
     N_components_(N_components),
     extended_(isExtendedHelper(pdf_.get(), extended))
{
   if (extended == Extended::Auto) {
      if (extended_) {
         oocoutI(nullptr, Minimization)
            << "in RooAbsL ctor: p.d.f. provides expected number of events, including extended term in likelihood."
            << std::endl;
      }
   }
}
 
/// Constructor that clones the pdf/data and owns those cloned copies.
///
/// This constructor is used for classes that need a pdf/data clone (RooBinnedL and RooUnbinnedL).
///
/// \param in Struct containing raw pointers to the pdf and dataset that are to be cloned.
/// \param N_events The number of events in this likelihood's dataset.
/// \param N_components The number of components in the likelihood.
/// \param extended Set extended term calculation on, off or use Extended::Auto to determine automatically based on the
/// pdf whether to activate or not.
RooAbsL::RooAbsL(RooAbsL::ClonePdfData in, std::size_t N_events, std::size_t N_components, Extended extended)
   : RooAbsL(in.ownedPdf ? std::move(in.ownedPdf)
                         : std::unique_ptr<RooAbsPdf>(static_cast<RooAbsPdf *>(in.pdf->cloneTree())),
             std::shared_ptr<RooAbsData>(static_cast<RooAbsData *>(in.data->Clone())), N_events, N_components, extended)
{
   initClones(*in.pdf, *in.data);
}
 
/// Constructor that does not clone pdf/data and uses the shared_ptr aliasing constructor to make it non-owning.
///
/// This constructor is used for classes where a reference to the external pdf/dataset is good enough (RooSumL and
/// RooSubsidiaryL).
///
/// \param inpdf Raw pointer to the pdf.
/// \param indata Raw pointer to the dataset.
/// \param N_events The number of events in this likelihood's dataset.
/// \param N_components The number of components in the likelihood.
/// \param extended Set extended term calculation on, off or use Extended::Auto to determine automatically based on the
/// pdf whether to activate or not.
RooAbsL::RooAbsL(RooAbsPdf *inpdf, RooAbsData *indata, std::size_t N_events, std::size_t N_components,
                 Extended extended)
   : RooAbsL({std::shared_ptr<RooAbsPdf>(nullptr), inpdf}, {std::shared_ptr<RooAbsData>(nullptr), indata}, N_events,
             N_components, extended)
{
}
 
RooAbsL::RooAbsL(const RooAbsL &other)
   : pdf_(other.pdf_),
     data_(other.data_),
     N_events_(other.N_events_),
     N_components_(other.N_components_),
     extended_(other.extended_),
     sim_count_(other.sim_count_)
{
   // it can never be one, since we just copied the shared_ptr; if it is, something really weird is going on; also they
   // must be equal (usually either zero or two)
   assert((pdf_.use_count() != 1) && (data_.use_count() != 1) && (pdf_.use_count() == data_.use_count()));
   if ((pdf_.use_count() > 1) && (data_.use_count() > 1)) {
      pdf_.reset(static_cast<RooAbsPdf *>(other.pdf_->cloneTree()));
      data_.reset(static_cast<RooAbsData *>(other.data_->Clone()));
      initClones(*other.pdf_, *other.data_);
   }
}
 
void RooAbsL::initClones(RooAbsPdf &inpdf, RooAbsData &indata)
{
   // ******************************************************************
   // *** PART 1 *** Clone incoming pdf, attach to each other *
   // ******************************************************************
 
   // Attach FUNC to data set
   auto _funcObsSet = pdf_->getObservables(indata);
 
   if (pdf_->getAttribute("BinnedLikelihood")) {
      pdf_->setAttribute("BinnedLikelihoodActive");
   }
 
   // Reattach FUNC to original parameters
   std::unique_ptr<RooArgSet> origParams{inpdf.getParameters(indata)};
   pdf_->recursiveRedirectServers(*origParams);
 
   // Store normalization set
   normSet_ = std::make_unique<RooArgSet>();
   indata.get()->snapshot(*normSet_, false);
 
   // Expand list of observables with any observables used in parameterized ranges
   for (const auto realDep : *_funcObsSet) {
      auto realDepRLV = dynamic_cast<RooAbsRealLValue *>(realDep);
      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
      auto realReal = dynamic_cast<RooRealVar *>(arg);
      if (!realReal) {
         continue;
      }
      auto 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)) {
         oocoutE(nullptr, InputArguments) << "RooAbsL: 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)) {
         oocoutE(nullptr, InputArguments)
            << "RooAbsL: ERROR maximum of FUNC observable " << arg->GetName() << " is larger than that of "
            << arg->GetName() << " in the dataset" << std::endl;
         RooErrorHandler::softAbort();
         return;
      }
   }
 
   // ******************************************************************
   // *** PART 3 *** Make adjustments for fit ranges, if specified     *
   // ******************************************************************
 
   // TODO
 
   // Jonas R.: The following code is commented out, because the functionality
   // to mask out-of-range entries with `RooDataHist::cacheValidEntries` has
   // been removed from the RooDataHist. If you want to implement ranged fits
   // properly, please create a RooDataHist for the requested range with
   // `RooDataHist::reduce`.
 
   //// If dataset is binned, activate caching of bins that are invalid because they're outside the
   //// updated range definition (WVE need to add virtual interface here)
   // RooDataHist *tmph = dynamic_cast<RooDataHist *>(data_.get());
   // if (tmph) {
   // tmph->cacheValidEntries();
   //}
 
   // This is deferred from part 2 - but must happen after part 3 - otherwise invalid bins cannot be properly marked in
   // cacheValidEntries
   data_->attachBuffers(*_funcObsSet);
 
   // *********************************************************************
   // *** PART 4 *** Adjust normalization range for projected observables *
   // *********************************************************************
 
   // TODO
 
   // *********************************************************************
   // *** PART 4 *** Finalization and activation of optimization          *
   // *********************************************************************
 
   // optimization steps (copied from ROATS::optimizeCaching)
 
   pdf_->getVal(normSet_.get());
   // Set value caching mode for all nodes that depend on any of the observables to ADirty
   pdf_->optimizeCacheMode(*_funcObsSet);
   // Disable propagation of dirty state flags for observables
   data_->setDirtyProp(false);
 
   // Disable reading of observables that are not used
   data_->optimizeReadingWithCaching(*pdf_, RooArgSet(), RooArgSet());
}
 
std::unique_ptr<RooArgSet> RooAbsL::getParameters()
{
   return std::unique_ptr<RooArgSet>{pdf_->getParameters(*data_)};
}
 
void RooAbsL::constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt)
{
   if (data_.get()->hasFilledCache() && opcode == RooAbsArg::Activate) {
      opcode = RooAbsArg::ValueChange;
   }
 
   switch (opcode) {
   case RooAbsArg::Activate:
      oocxcoutI((TObject *)nullptr, Optimization)
         << "RooAbsL::constOptimizeTestStatistic(" << GetName()
         << ") optimizing evaluation of test statistic by finding all nodes in p.d.f that depend exclusively"
         << " on observables and constant parameters and precalculating their values" << std::endl;
      ConstantTermsOptimizer::enableConstantTermsOptimization(pdf_.get(), normSet_.get(), data_.get(),
                                                              doAlsoTrackingOpt);
      break;
 
   case RooAbsArg::DeActivate:
      oocxcoutI((TObject *)nullptr, Optimization)
         << "RooAbsL::constOptimizeTestStatistic(" << GetName()
         << ") deactivating optimization of constant terms in test statistic" << std::endl;
      ConstantTermsOptimizer::disableConstantTermsOptimization(pdf_.get(), normSet_.get(), data_.get());
      break;
 
   case RooAbsArg::ConfigChange:
      oocxcoutI((TObject *)nullptr, Optimization)
         << "RooAbsL::constOptimizeTestStatistic(" << GetName()
         << ") one ore more parameter were changed from constant to floating or vice versa, "
         << "re-evaluating constant term optimization" << std::endl;
      ConstantTermsOptimizer::disableConstantTermsOptimization(pdf_.get(), normSet_.get(), data_.get());
      ConstantTermsOptimizer::enableConstantTermsOptimization(pdf_.get(), normSet_.get(), data_.get(),
                                                              doAlsoTrackingOpt);
      break;
 
   case RooAbsArg::ValueChange:
      oocxcoutI((TObject *)nullptr, Optimization)
         << "RooAbsL::constOptimizeTestStatistic(" << GetName()
         << ") the value of one ore more constant parameter were changed re-evaluating constant term optimization"
         << std::endl;
      // Request a forcible cache update of all cached nodes
      data_.get()->store()->forceCacheUpdate();
      break;
   }
}
 
std::string RooAbsL::GetName() const
{
   std::string output("likelihood of pdf ");
   output.append(pdf_->GetName());
   return output;
}
 
std::string RooAbsL::GetTitle() const
{
   std::string output("likelihood of pdf ");
   output.append(pdf_->GetTitle());
   return output;
}
 
std::size_t RooAbsL::numDataEntries() const
{
   return static_cast<std::size_t>(data_->numEntries());
}
 
} // namespace TestStatistics
} // namespace RooFit