diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index 089849d4de68d..b61c89b77d5f0 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -233,6 +233,19 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore RooFitLegacy/RooNameSet.h RooFitLegacy/RooSetPair.h RooFitLegacy/RooTreeData.h + TestStatistics/ConstantTermsOptimizer.h + TestStatistics/LikelihoodGradientWrapper.h + TestStatistics/LikelihoodWrapper.h + TestStatistics/LikelihoodSerial.h + TestStatistics/MinuitFcnGrad.h + TestStatistics/RooAbsL.h + TestStatistics/RooBinnedL.h + TestStatistics/RooSubsidiaryL.h + TestStatistics/RooSumL.h + TestStatistics/RooRealL.h + TestStatistics/RooUnbinnedL.h + TestStatistics/optional_parameter_types.h + TestStatistics/buildLikelihood.h SOURCES src/BidirMMapPipe.cxx src/BidirMMapPipe.h @@ -449,6 +462,19 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore src/RooFitLegacy/RooMultiCatIter.cxx src/RooFitLegacy/RooNameSet.cxx src/RooFitLegacy/RooSetPair.cxx + src/TestStatistics/ConstantTermsOptimizer.cxx + src/TestStatistics/LikelihoodGradientWrapper.cxx + src/TestStatistics/LikelihoodWrapper.cxx + src/TestStatistics/LikelihoodSerial.cxx + src/TestStatistics/MinuitFcnGrad.cxx + src/TestStatistics/RooAbsL.cxx + src/TestStatistics/RooBinnedL.cxx + src/TestStatistics/RooSubsidiaryL.cxx + src/TestStatistics/RooSumL.cxx + src/TestStatistics/RooRealL.cxx + src/TestStatistics/RooUnbinnedL.cxx + src/TestStatistics/optional_parameter_types.cxx + src/TestStatistics/buildLikelihood.cxx DICTIONARY_OPTIONS "-writeEmptyRootPCM" LIBRARIES diff --git a/roofit/roofitcore/inc/LinkDef.h b/roofit/roofitcore/inc/LinkDef.h index 9b54f2b7e5ee1..2c2d14b7770e7 100644 --- a/roofit/roofitcore/inc/LinkDef.h +++ b/roofit/roofitcore/inc/LinkDef.h @@ -338,6 +338,7 @@ #pragma link C++ class RooUnitTest+ ; #pragma link C++ class RooMinimizer+ ; #pragma link C++ class RooMinimizerFcn+ ; +#pragma link C++ class RooFit::TestStatistics::RooRealL+ ; #pragma link C++ class RooAbsMoment+ ; #pragma link C++ class RooMoment+ ; #pragma link C++ class RooFirstMoment+ ; diff --git a/roofit/roofitcore/inc/RooAbsData.h b/roofit/roofitcore/inc/RooAbsData.h index 7a5b19a980fa4..6657583f9c1ad 100644 --- a/roofit/roofitcore/inc/RooAbsData.h +++ b/roofit/roofitcore/inc/RooAbsData.h @@ -47,6 +47,12 @@ class RooFormulaVar; namespace RooBatchCompute{ struct RunContext; } +namespace RooFit { +namespace TestStatistics { +class RooAbsL; +struct ConstantTermsOptimizer; +} +} // Writes a templated constructor for compatibility with ROOT builds using the @@ -321,6 +327,9 @@ class RooAbsData : public TNamed, public RooPrintable { friend class RooAbsReal ; friend class RooAbsOptTestStatistic ; friend class RooAbsCachedPdf ; + friend struct RooFit::TestStatistics::ConstantTermsOptimizer; + // for access into copied dataset: + friend class RooFit::TestStatistics::RooAbsL; virtual void cacheArgs(const RooAbsArg* owner, RooArgSet& varSet, const RooArgSet* nset=0, Bool_t skipZeroWeights=kFALSE) ; virtual void resetCache() ; diff --git a/roofit/roofitcore/inc/RooDataHist.h b/roofit/roofitcore/inc/RooDataHist.h index 59d4652932ea7..abcd1130f160a 100644 --- a/roofit/roofitcore/inc/RooDataHist.h +++ b/roofit/roofitcore/inc/RooDataHist.h @@ -35,6 +35,11 @@ class RooAbsArg; class RooCategory ; class RooPlot; class RooAbsLValue ; +namespace RooFit { +namespace TestStatistics { +class RooAbsL; +} +} class RooDataHist : public RooAbsData, public RooDirItem { public: @@ -222,6 +227,8 @@ class RooDataHist : public RooAbsData, public RooDirItem { friend class RooAbsCachedPdf ; friend class RooAbsCachedReal ; friend class RooDataHistSliceIter ; + // for access into copied dataset: + friend class RooFit::TestStatistics::RooAbsL; std::size_t calcTreeIndex(const RooAbsCollection& coords, bool fast) const; /// Legacy overload to calculate the tree index from the current value of `_vars`. diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h index 835e71b0aa791..6d590a983d547 100644 --- a/roofit/roofitcore/inc/RooMinimizer.h +++ b/roofit/roofitcore/inc/RooMinimizer.h @@ -29,8 +29,8 @@ #include "RooArgList.h" // cannot just use forward decl due to default argument in lastMinuitFit -#include "RooMinimizerFcn.h" -#include "RooGradMinimizerFcn.h" +#include "TestStatistics/RooAbsL.h" +#include "TestStatistics/MinuitFcnGrad.h" #include "RooSentinel.h" #include "RooMsgService.h" @@ -45,13 +45,23 @@ class RooRealVar ; class RooArgSet ; class TH2F ; class RooPlot ; +namespace RooFit { +namespace TestStatistics { +class LikelihoodSerial; +class LikelihoodGradientSerial; +} +} // namespace RooFit class RooMinimizer : public TObject { public: - enum class FcnMode { classic, gradient }; + enum class FcnMode { classic, gradient, generic_wrapper }; explicit RooMinimizer(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic); static std::unique_ptr create(RooAbsReal &function, FcnMode fcnMode = FcnMode::classic); + template + static std::unique_ptr create(std::shared_ptr likelihood); + ~RooMinimizer() override; enum Strategy { Speed=0, Balance=1, Robustness=2 } ; @@ -126,6 +136,11 @@ class RooMinimizer : public TObject { bool fitFcn() const; private: + template + RooMinimizer(std::shared_ptr likelihood, + LikelihoodWrapperT* /* used only for template deduction */ = static_cast(nullptr), + LikelihoodGradientWrapperT* /* used only for template deduction */ = static_cast(nullptr)); + Int_t _printLevel = 1; Int_t _status = -99; Bool_t _profile = kFALSE; @@ -150,4 +165,48 @@ class RooMinimizer : public TObject { ClassDefOverride(RooMinimizer,0) // RooFit interface to ROOT::Fit::Fitter } ; + +template +RooMinimizer::RooMinimizer(std::shared_ptr likelihood, LikelihoodWrapperT* /* value unused */, + LikelihoodGradientWrapperT* /* value unused */) : + _fcnMode(FcnMode::generic_wrapper) +{ + RooSentinel::activate(); + + if (_theFitter) + delete _theFitter; + _theFitter = new ROOT::Fit::Fitter; + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + setEps(1.0); // default tolerance + + _fcn = RooFit::TestStatistics::MinuitFcnGrad::create(likelihood, this, _theFitter->Config().ParamsSettings(), _verbose); + + // default max number of calls + _theFitter->Config().MinimizerOptions().SetMaxIterations(500 * _fcn->getNDim()); + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500 * _fcn->getNDim()); + + // Shut up for now + setPrintLevel(-1); + + // Use +0.5 for 1-sigma errors + setErrorLevel(likelihood->defaultErrorLevel()); + + // Declare our parameters to MINUIT + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), _fcn->getOptConst(), _verbose); + + // Now set default verbosity + if (RooMsgService::instance().silentMode()) { + setPrintLevel(-1); + } else { + setPrintLevel(1); + } +} + +// static function +template +std::unique_ptr RooMinimizer::create(std::shared_ptr likelihood) { + return std::make_unique(likelihood, static_cast(nullptr), + static_cast(nullptr)); +} + #endif diff --git a/roofit/roofitcore/inc/RooNLLVar.h b/roofit/roofitcore/inc/RooNLLVar.h index daecdc9b4a8df..78019a9066a40 100644 --- a/roofit/roofitcore/inc/RooNLLVar.h +++ b/roofit/roofitcore/inc/RooNLLVar.h @@ -62,6 +62,14 @@ class RooNLLVar : public RooAbsOptTestStatistic { using ComputeResult = std::pair, double>; + static RooNLLVar::ComputeResult computeBatchedFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, + std::unique_ptr &evalData, + RooArgSet *normSet, bool weightSq, std::size_t stepSize, + std::size_t firstEvent, std::size_t lastEvent); + static RooNLLVar::ComputeResult computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, RooArgSet *normSet, + bool weightSq, std::size_t stepSize, std::size_t firstEvent, + std::size_t lastEvent); + protected: virtual Bool_t processEmptyDataSets() const { return _extended ; } diff --git a/roofit/roofitcore/inc/TestStatistics/ConstantTermsOptimizer.h b/roofit/roofitcore/inc/TestStatistics/ConstantTermsOptimizer.h new file mode 100644 index 0000000000000..ff303c07d50d1 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/ConstantTermsOptimizer.h @@ -0,0 +1,40 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_ConstantTermsOptimizer +#define ROOT_ROOFIT_TESTSTATISTICS_ConstantTermsOptimizer + +// forward declarations +class RooAbsReal; +class RooArgSet; +class RooAbsData; + +namespace RooFit { +namespace TestStatistics { + +// this is a class only for convenience: it saves multiple friend definitions in RooAbsData for otherwise free functions +struct ConstantTermsOptimizer { + static void enableConstantTermsOptimization(RooAbsReal *function, RooArgSet *norm_set, RooAbsData *dataset, + bool applyTrackingOpt); + static void optimizeCaching(RooAbsReal *function, RooArgSet *norm_set, RooArgSet* observables, RooAbsData *dataset); + static void disableConstantTermsOptimization(RooAbsReal *function, RooArgSet *norm_set, RooArgSet* observables, RooAbsData *dataset); + static RooArgSet requiredExtraObservables(); +}; + +} +} + +#endif // ROOT_ROOFIT_TESTSTATISTICS_ConstantTermsOptimizer diff --git a/roofit/roofitcore/inc/TestStatistics/LikelihoodGradientWrapper.h b/roofit/roofitcore/inc/TestStatistics/LikelihoodGradientWrapper.h new file mode 100644 index 0000000000000..44bb4397945ad --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/LikelihoodGradientWrapper.h @@ -0,0 +1,73 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_LikelihoodGradientWrapper +#define ROOT_ROOFIT_TESTSTATISTICS_LikelihoodGradientWrapper + +#include +#include +#include "Math/MinimizerOptions.h" + +#include +#include // shared_ptr + +// forward declaration +class RooMinimizer; + +namespace RooFit { +namespace TestStatistics { + +// forward declaration +class RooAbsL; +struct WrapperCalculationCleanFlags; + +class LikelihoodGradientWrapper { +public: + LikelihoodGradientWrapper(std::shared_ptr likelihood, std::shared_ptr calculation_is_clean, std::size_t N_dim, RooMinimizer* minimizer); + virtual ~LikelihoodGradientWrapper() = default; + virtual LikelihoodGradientWrapper* clone() const = 0; + + virtual void fillGradient(double *grad) = 0; + + /// Synchronize minimizer settings with calculators in child classes. + virtual void synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions &options); + virtual void synchronizeParameterSettings(const std::vector ¶meter_settings); + virtual void synchronizeParameterSettings(ROOT::Math::IMultiGenFunction* function, const std::vector ¶meter_settings) = 0; + /// Minuit passes in parameter values that may not conform to RooFit internal standards (like applying range clipping), + /// but that the specific calculator does need. This function can be implemented to receive these Minuit-internal values. + virtual void updateMinuitInternalParameterValues(const std::vector& minuit_internal_x); + virtual void updateMinuitExternalParameterValues(const std::vector& minuit_external_x); + + /// \brief Implement usesMinuitInternalValues to return true when you want Minuit to send this class Minuit-internal values, or return false when you want "regular" Minuit-external values. + /// + /// Minuit internally uses a transformed parameter space to graciously handle externally mandated parameter range + /// boundaries. Transformation from Minuit-internal to external (i.e. "regular") parameters is done using + /// trigonometric functions that in some cases can cause a few bits of precision loss with respect to the original + /// parameter values. To circumvent this, Minuit also allows external gradient providers (like LikelihoodGradientWrapper) + /// to take the Minuit-internal parameter values directly, without transformation. This way, the gradient provider + /// (e.g. the implementation of this class) can handle transformation manually, possibly with higher precision. + virtual bool usesMinuitInternalValues() = 0; + +protected: + std::shared_ptr likelihood_; + RooMinimizer * minimizer_; + std::shared_ptr calculation_is_clean_; +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_LikelihoodGradientWrapper diff --git a/roofit/roofitcore/inc/TestStatistics/LikelihoodSerial.h b/roofit/roofitcore/inc/TestStatistics/LikelihoodSerial.h new file mode 100644 index 0000000000000..720483e3bfc46 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/LikelihoodSerial.h @@ -0,0 +1,52 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_LikelihoodSerial +#define ROOT_ROOFIT_LikelihoodSerial + +#include +#include "RooArgList.h" + +#include "Math/MinimizerOptions.h" + +#include + +namespace RooFit { +namespace TestStatistics { + +class LikelihoodSerial : public LikelihoodWrapper { +public: + LikelihoodSerial(std::shared_ptr likelihood, std::shared_ptr calculation_is_clean); + inline LikelihoodSerial *clone() const override { return new LikelihoodSerial(*this); } + + void initVars(); + + void evaluate() override; + inline ROOT::Math::KahanSum getResult() const override { return result; } + +private: + ROOT::Math::KahanSum result; + + RooArgList _vars; // Variables + RooArgList _saveVars; // Copy of variables + + LikelihoodType likelihood_type; +}; + +} +} + +#endif // ROOT_ROOFIT_LikelihoodSerial diff --git a/roofit/roofitcore/inc/TestStatistics/LikelihoodWrapper.h b/roofit/roofitcore/inc/TestStatistics/LikelihoodWrapper.h new file mode 100644 index 0000000000000..45973089485e3 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/LikelihoodWrapper.h @@ -0,0 +1,107 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_LikelihoodWrapper +#define ROOT_ROOFIT_TESTSTATISTICS_LikelihoodWrapper + +#include "RooArgSet.h" +#include "RooAbsArg.h" // enum ConstOpCode + +#include +#include "Math/MinimizerOptions.h" + +#include // shared_ptr +#include + +// forward declaration +class RooMinimizer; + +namespace RooFit { +namespace TestStatistics { + +// forward declaration +class RooAbsL; +struct WrapperCalculationCleanFlags; + +enum class LikelihoodType { + unbinned, + binned, + subsidiary, + sum +}; + +/// Previously, offsetting was only implemented for RooNLLVar components of a likelihood, +/// not for RooConstraintSum terms. To emulate this behavior, use OffsettingMode::legacy. To +/// also offset the RooSubsidiaryL component (equivalent of RooConstraintSum) of RooSumL +/// likelihoods, use OffsettingMode::full. +enum class OffsettingMode { + legacy, + full +}; + +class LikelihoodWrapper { +public: + LikelihoodWrapper(std::shared_ptr likelihood, std::shared_ptr calculation_is_clean); + virtual ~LikelihoodWrapper() = default; + virtual LikelihoodWrapper* clone() const = 0; + + /// \brief Triggers (possibly asynchronous) evaluation of the likelihood + /// + /// In parallel strategies, it may be advantageous to allow a calling process to continue on with other tasks while + /// the calculation is offloaded to another process or device, like a GPU. For this reason, evaluate() does not + /// return the result, this is done in getResult(). + virtual void evaluate() = 0; + /// \brief Return the latest result of a likelihood evaluation. + /// + /// Returns the result that was stored after calling evaluate(). It is up to the implementer to make sure the stored + /// value represents the most recent evaluation call, e.g. by using a mutex. + virtual ROOT::Math::KahanSum getResult() const = 0; + + /// Synchronize minimizer settings with calculators in child classes + virtual void synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions & options); + virtual void synchronizeParameterSettings(const std::vector ¶meter_settings); + /// Minuit passes in parameter values that may not conform to RooFit internal standards (like applying range clipping), + /// but that the specific calculator does need. This function can be implemented to receive these Minuit-internal values: + virtual void updateMinuitInternalParameterValues(const std::vector& minuit_internal_x); + virtual void updateMinuitExternalParameterValues(const std::vector& minuit_external_x); + + // The following functions are necessary from MinuitFcnGrad to reach likelihood properties: + void constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt); + double defaultErrorLevel() const; + virtual std::string GetName() const; + virtual std::string GetTitle() const; + inline virtual bool isOffsetting() const { return do_offset_; } + virtual void enableOffsetting(bool flag); + void setOffsettingMode(OffsettingMode mode); + inline ROOT::Math::KahanSum offset() const { return offset_; } + void setApplyWeightSquared(bool flag); + +protected: + std::shared_ptr likelihood_; + std::shared_ptr calculation_is_clean_; + + bool do_offset_ = false; + ROOT::Math::KahanSum offset_; + ROOT::Math::KahanSum offset_save_ = 0; //! + OffsettingMode offsetting_mode_ = OffsettingMode::legacy; + ROOT::Math::KahanSum applyOffsetting(ROOT::Math::KahanSum current_value); + void swapOffsets(); +}; + +} +} + +#endif // ROOT_ROOFIT_TESTSTATISTICS_LikelihoodWrapper diff --git a/roofit/roofitcore/inc/TestStatistics/MinuitFcnGrad.h b/roofit/roofitcore/inc/TestStatistics/MinuitFcnGrad.h new file mode 100644 index 0000000000000..5dca9126adc2a --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/MinuitFcnGrad.h @@ -0,0 +1,166 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_MinuitFcnGrad +#define ROOT_ROOFIT_TESTSTATISTICS_MinuitFcnGrad + +#include "RooArgList.h" +#include "RooRealVar.h" +#include "TestStatistics/RooAbsL.h" +#include "TestStatistics/LikelihoodWrapper.h" +#include "TestStatistics/LikelihoodGradientWrapper.h" +#include "RooAbsMinimizerFcn.h" + +#include +#include +#include "Math/IFunction.h" // ROOT::Math::IMultiGradFunction + +// forward declaration +class RooAbsReal; +class RooMinimizer; + +namespace RooFit { +namespace TestStatistics { + +// forward declaration +class LikelihoodSerial; +class LikelihoodGradientSerial; + +/// For communication with wrappers, an instance of this struct must be shared between them and MinuitFcnGrad. It keeps +/// track of what has been evaluated for the current parameter set provided by Minuit. +struct WrapperCalculationCleanFlags { + // indicate whether that part has been calculated since the last parameter update + bool likelihood = false; + bool gradient = false; + + void set_all(bool value) + { + likelihood = value; + gradient = value; + } +}; + +class MinuitFcnGrad : public ROOT::Math::IMultiGradFunction, public RooAbsMinimizerFcn { +public: + // factory + template + static MinuitFcnGrad * + create(const std::shared_ptr &likelihood, RooMinimizer *context, + std::vector ¶meters, bool verbose = false); + + inline ROOT::Math::IMultiGradFunction *Clone() const override { return new MinuitFcnGrad(*this); } + + /// Overridden from RooAbsMinimizerFcn to include gradient strategy synchronization. + Bool_t Synchronize(std::vector ¶meter_settings, Bool_t optConst, + Bool_t verbose = kFALSE) override; + + // used inside Minuit: + inline bool returnsInMinuit2ParameterSpace() const override { return gradient->usesMinuitInternalValues(); } + + inline void setOptimizeConstOnFunction(RooAbsArg::ConstOpCode opcode, Bool_t doAlsoTrackingOpt) override + { + likelihood->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); + } + +private: + /// IMultiGradFunction override necessary for Minuit + double DoEval(const double *x) const override; + +public: + /// IMultiGradFunction override necessary for Minuit + void Gradient(const double *x, double *grad) const override; + + /// Part of IMultiGradFunction interface, used widely both in Minuit and in RooFit. + inline unsigned int NDim() const override { return _nDim; } + + inline std::string getFunctionName() const override { return likelihood->GetName(); } + + inline std::string getFunctionTitle() const override { return likelihood->GetTitle(); } + + inline void setOffsetting(Bool_t flag) override { likelihood->enableOffsetting(flag); } + +private: + template + MinuitFcnGrad( + const std::shared_ptr &_likelihood, RooMinimizer *context, + std::vector ¶meters, bool verbose, + LikelihoodWrapperT * /* used only for template deduction */ = static_cast(nullptr), + LikelihoodGradientWrapperT * /* used only for template deduction */ = + static_cast(nullptr)); + + /// This override should not be used in this class, so it throws. + double DoDerivative(const double *x, unsigned int icoord) const override; + + bool syncParameterValuesFromMinuitCalls(const double *x, bool minuit_internal) const; + + // members + std::shared_ptr likelihood; + std::shared_ptr gradient; + +public: + mutable std::shared_ptr calculation_is_clean; + +private: + mutable std::vector minuit_internal_x_; + mutable std::vector minuit_external_x_; + +public: + mutable bool minuit_internal_roofit_x_mismatch_ = false; +}; + + +/// \param[in] context RooMinimizer that creates and owns this class. +/// \param[in] parameters The vector of ParameterSettings objects that describe the parameters used in the Minuit +/// Fitter. Note that these must match the set used in the Fitter used by \p context! It can be passed in from +/// RooMinimizer with fitter()->Config().ParamsSettings(). +template +MinuitFcnGrad::MinuitFcnGrad(const std::shared_ptr &_likelihood, RooMinimizer *context, + std::vector ¶meters, bool verbose, + LikelihoodWrapperT * /* value unused */, LikelihoodGradientWrapperT * /* value unused */) + : RooAbsMinimizerFcn(RooArgList(*_likelihood->getParameters()), context, verbose), minuit_internal_x_(NDim(), 0), + minuit_external_x_(NDim(), 0) +{ + synchronizeParameterSettings(parameters, kTRUE, verbose); + + calculation_is_clean = std::make_shared(); + likelihood = std::make_shared(_likelihood, calculation_is_clean /*, _context*/); + gradient = std::make_shared(_likelihood, calculation_is_clean, getNDim(), _context); + + likelihood->synchronizeParameterSettings(parameters); + gradient->synchronizeParameterSettings(this, parameters); + + // Note: can be different than RooGradMinimizerFcn, where default options are passed + // (ROOT::Math::MinimizerOptions::DefaultStrategy() and ROOT::Math::MinimizerOptions::DefaultErrorDef()) + likelihood->synchronizeWithMinimizer(ROOT::Math::MinimizerOptions()); + gradient->synchronizeWithMinimizer(ROOT::Math::MinimizerOptions()); +} + +// static function +template +MinuitFcnGrad * +MinuitFcnGrad::create(const std::shared_ptr &likelihood, RooMinimizer *context, + std::vector ¶meters, bool verbose) +{ + return new MinuitFcnGrad(likelihood, context, parameters, verbose, static_cast(nullptr), + static_cast(nullptr)); +} + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_MinuitFcnGrad diff --git a/roofit/roofitcore/inc/TestStatistics/RooAbsL.h b/roofit/roofitcore/inc/TestStatistics/RooAbsL.h new file mode 100644 index 0000000000000..68881b6065e48 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooAbsL.h @@ -0,0 +1,148 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooAbsL +#define ROOT_ROOFIT_TESTSTATISTICS_RooAbsL + +#include "RooArgSet.h" +#include "RooAbsArg.h" // enum ConstOpCode + +#include "Math/Util.h" // KahanSum + +#include // std::size_t +#include +#include + +// forward declarations +class RooAbsPdf; +class RooAbsData; + +namespace RooFit { +namespace TestStatistics { + +class RooAbsL { +public: + enum class Extended { Auto, Yes, No }; + static bool isExtendedHelper(RooAbsPdf *pdf, Extended extended); + + /// Convenience wrapper class used to distinguish between pdf/data owning and non-owning constructors. + struct ClonePdfData { + RooAbsPdf *pdf; + RooAbsData *data; + }; + +private: + RooAbsL(std::shared_ptr pdf, std::shared_ptr data, std::size_t N_events, + std::size_t N_components, Extended extended); + +public: + RooAbsL(RooAbsPdf *pdf, RooAbsData *data, std::size_t N_events, std::size_t N_components, + Extended extended = Extended::Auto); + RooAbsL(ClonePdfData in, std::size_t N_events, std::size_t N_components, Extended extended = Extended::Auto); + RooAbsL(const RooAbsL &other); + virtual ~RooAbsL() = default; + + void initClones(RooAbsPdf &inpdf, RooAbsData &indata); + + /// A part of some range delimited by two fractional points between 0 and 1 (inclusive). + struct Section { + Section(double begin, double end) : begin_fraction(begin), end_fraction(end) + { + if ((begin > end) || (begin < 0) || (end > 1)) { + throw std::domain_error("Invalid input values for section; begin must be >= 0, end <= 1 and begin < end."); + } + } + + Section(const Section §ion) = default; + + std::size_t begin(std::size_t N_total) const { return static_cast(N_total * begin_fraction); } + + std::size_t end(std::size_t N_total) const + { + if (end_fraction == 1) { + return N_total; + } else { + return static_cast(N_total * end_fraction); + } + } + + double begin_fraction; + double end_fraction; + }; + + /* + * \brief Evaluate (part of) the likelihood over a given range of events and components + * + * A fractional event range is used because components may have different numbers of events. For a + * multi-component RooSumL, for instance, this means the caller need not indicate for each component which event + * ranges they want to evaluate, but can just pass one overall fractional range. + * + * \param[in] events The fractional event range. + * \param[in] components_begin The first component to be calculated. + * \param[in] components_end The *exclusive* upper limit to the range of components to be calculated, i.e. the + * component *before this one* is the last to be included. \return The value of part of the negative log likelihood, + * returned as a KahanSum object which also includes a carry term. + */ + virtual ROOT::Math::KahanSum + evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) = 0; + + // necessary from MinuitFcnGrad to reach likelihood properties: + virtual RooArgSet *getParameters(); + + /// \brief Interface function signaling a request to perform constant term optimization. + /// + /// The default implementation takes no action other than to forward the calls to all servers. May be overridden in + /// likelihood classes without a cached dataset, like RooSubsidiaryL. + virtual void constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt); + + virtual std::string GetName() const; + virtual std::string GetTitle() const; + + // necessary in RooMinimizer (via LikelihoodWrapper) + inline virtual double defaultErrorLevel() const { return 0.5; } + + // necessary in LikelihoodJob + /// Number of dataset entries. Typically equal to the number of dataset events, except in RooSubsidiaryL, which has + /// no events. + virtual std::size_t numDataEntries() const; + inline std::size_t getNEvents() const { return N_events_; } + inline std::size_t getNComponents() const { return N_components_; } + inline bool isExtended() const { return extended_; } + inline void setSimCount(std::size_t value) { sim_count_ = value; } + +protected: + // Note: pdf_ and data_ can be constructed in two ways, one of which implies ownership and the other does not. + // Inspired by this: https://stackoverflow.com/a/61227865/1199693. + // The owning variant is used for classes that need a pdf/data clone (RooBinnedL and RooUnbinnedL), whereas the + // non-owning version is used for when a reference to the external pdf/dataset is good enough (RooSumL). + // This means that pdf_ and data_ are not meant to actually be shared! If there were a unique_ptr with optional + // ownership, we would have used that instead. + std::shared_ptr pdf_; + std::shared_ptr data_; + std::unique_ptr normSet_; // Pointer to set with observables used for normalization + + std::size_t N_events_ = 1; + std::size_t N_components_ = 1; + + bool extended_ = false; + + std::size_t sim_count_ = 1; // Total number of component p.d.f.s in RooSimultaneous (if any) +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooAbsL diff --git a/roofit/roofitcore/inc/TestStatistics/RooBinnedL.h b/roofit/roofitcore/inc/TestStatistics/RooBinnedL.h new file mode 100644 index 0000000000000..2a72a5cccfe89 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooBinnedL.h @@ -0,0 +1,45 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooBinnedL +#define ROOT_ROOFIT_TESTSTATISTICS_RooBinnedL + +#include +#include "RooAbsReal.h" +#include + +// forward declarations +class RooAbsPdf; +class RooAbsData; + +namespace RooFit { +namespace TestStatistics { + +class RooBinnedL : public RooAbsL { +public: + RooBinnedL(RooAbsPdf *pdf, RooAbsData *data); + ROOT::Math::KahanSum + evaluatePartition(Section bins, std::size_t components_begin, std::size_t components_end) override; + +private: + mutable bool _first = true; //! + mutable std::vector _binw; //! +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooBinnedL diff --git a/roofit/roofitcore/inc/TestStatistics/RooRealL.h b/roofit/roofitcore/inc/TestStatistics/RooRealL.h new file mode 100644 index 0000000000000..3effc137cd80a --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooRealL.h @@ -0,0 +1,60 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooRealL +#define ROOT_ROOFIT_TESTSTATISTICS_RooRealL + +#include "RooAbsReal.h" +#include "RooSetProxy.h" + +#include "Rtypes.h" // ClassDef, ClassImp + +#include // shared_ptr + +namespace RooFit { +namespace TestStatistics { + +class RooAbsL; + +class RooRealL : public RooAbsReal { +public: + RooRealL(const char *name, const char *title, std::shared_ptr likelihood); + RooRealL(const RooRealL &other, const char *name = 0); + + Double_t evaluate() const override; + inline TObject *clone(const char *newname) const override { return new RooRealL(*this, newname); } + + inline double globalNormalization() const + { + // Default value of global normalization factor is 1.0 + return 1.0; + } + + inline double getCarry() const { return eval_carry; } + inline Double_t defaultErrorLevel() const override { return 0.5; } + +private: + std::shared_ptr likelihood_; + mutable double eval_carry = 0; + RooSetProxy vars_proxy_; // sets up client-server connections + + ClassDefOverride(RooRealL, 0); +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooRealL diff --git a/roofit/roofitcore/inc/TestStatistics/RooSubsidiaryL.h b/roofit/roofitcore/inc/TestStatistics/RooSubsidiaryL.h new file mode 100644 index 0000000000000..d5a454df78319 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooSubsidiaryL.h @@ -0,0 +1,61 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooSubsidiaryL +#define ROOT_ROOFIT_TESTSTATISTICS_RooSubsidiaryL + +#include "TestStatistics/RooAbsL.h" +#include "RooArgList.h" +#include "RooArgSet.h" + +#include "Math/Util.h" // KahanSum + +namespace RooFit { +namespace TestStatistics { + +class RooSubsidiaryL : public RooAbsL { +public: + RooSubsidiaryL(const std::string &parent_pdf_name, const RooArgSet &pdfs, const RooArgSet ¶meter_set); + + ROOT::Math::KahanSum + evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) override; + inline RooArgSet *getParameters() override { return ¶meter_set_; } + inline std::string GetName() const override { return std::string("subsidiary_pdf_of_") + parent_pdf_name_; } + + inline std::string GetTitle() const override + { + return std::string("Subsidiary PDF set of simultaneous PDF ") + parent_pdf_name_; + } + + inline std::size_t numDataEntries() const override + { + // function only used in LikelihoodJob::evaluate, but this class must always be evaluated over Section(0,1), so + // not useful + return 0; + } + + void constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt) override; + +private: + std::string parent_pdf_name_; + RooArgList subsidiary_pdfs_{"subsidiary_pdfs"}; // Set of subsidiary PDF or "constraint" terms + RooArgSet parameter_set_{"parameter_set"}; // Set of parameters to which constraints apply +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooSubsidiaryL diff --git a/roofit/roofitcore/inc/TestStatistics/RooSumL.h b/roofit/roofitcore/inc/TestStatistics/RooSumL.h new file mode 100644 index 0000000000000..d8b30389f4700 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooSumL.h @@ -0,0 +1,52 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooSumL +#define ROOT_ROOFIT_TESTSTATISTICS_RooSumL + +#include "TestStatistics/RooAbsL.h" +#include "TestStatistics/optional_parameter_types.h" + +#include "Math/Util.h" // KahanSum + +#include + +namespace RooFit { +namespace TestStatistics { + +class RooSumL : public RooAbsL { +public: + RooSumL(RooAbsPdf *pdf, RooAbsData *data, std::vector> components, + RooAbsL::Extended extended = RooAbsL::Extended::Auto); + // Note: when above ctor is called without std::moving components, you get a really obscure error. Pass as + // std::move(components)! + + ROOT::Math::KahanSum + evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) override; + + // necessary only for legacy offsetting mode in LikelihoodWrapper; TODO: remove this if legacy mode is ever removed + ROOT::Math::KahanSum getSubsidiaryValue(); + + void constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt) override; + +private: + std::vector> components_; +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooSumL diff --git a/roofit/roofitcore/inc/TestStatistics/RooUnbinnedL.h b/roofit/roofitcore/inc/TestStatistics/RooUnbinnedL.h new file mode 100644 index 0000000000000..ae8cabfa9c550 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/RooUnbinnedL.h @@ -0,0 +1,54 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_RooUnbinnedL +#define ROOT_ROOFIT_TESTSTATISTICS_RooUnbinnedL + +#include + +// forward declarations +class RooAbsPdf; +class RooAbsData; +class RooArgSet; +namespace RooBatchCompute { +struct RunContext; +} + +namespace RooFit { +namespace TestStatistics { + +class RooUnbinnedL : public RooAbsL { +public: + RooUnbinnedL(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended = RooAbsL::Extended::Auto, + bool useBatchedEvaluations = false); + RooUnbinnedL(const RooUnbinnedL &other); + bool setApplyWeightSquared(bool flag); + + ROOT::Math::KahanSum + evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) override; + + void setUseBatchedEvaluations(bool flag); +private: + bool apply_weight_squared = false; // Apply weights squared? + mutable bool _first = true; //! + bool useBatchedEvaluations_ = false; + mutable std::unique_ptr evalData_; //! Struct to store function evaluation workspaces. +}; + +} // namespace TestStatistics +} // namespace RooFit + +#endif // ROOT_ROOFIT_TESTSTATISTICS_RooUnbinnedL diff --git a/roofit/roofitcore/inc/TestStatistics/buildLikelihood.h b/roofit/roofitcore/inc/TestStatistics/buildLikelihood.h new file mode 100644 index 0000000000000..d68ba253979b7 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/buildLikelihood.h @@ -0,0 +1,49 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_likelihood_builders +#define ROOT_ROOFIT_TESTSTATISTICS_likelihood_builders + +#include "TestStatistics/RooAbsL.h" +#include "TestStatistics/optional_parameter_types.h" + +#include + +// forward declarations +class RooAbsPdf; +class RooAbsData; + +namespace RooFit { +namespace TestStatistics { + +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended = RooAbsL::Extended::Auto, + ConstrainedParameters constrained_parameters = {}, ExternalConstraints external_constraints = {}, + GlobalObservables global_observables = {}, std::string global_observables_tag = {}); + +// delegating builder calls, for more convenient "optional" parameter passing +std::shared_ptr +buildLikelihood(RooAbsPdf* pdf, RooAbsData* data, ConstrainedParameters constrained_parameters); +std::shared_ptr buildLikelihood(RooAbsPdf* pdf, RooAbsData* data, ExternalConstraints external_constraints); +std::shared_ptr buildLikelihood(RooAbsPdf* pdf, RooAbsData* data, GlobalObservables global_observables); +std::shared_ptr buildLikelihood(RooAbsPdf* pdf, RooAbsData* data, std::string global_observables_tag); +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, ConstrainedParameters constrained_parameters, GlobalObservables global_observables); + +} +} + +#endif // ROOT_ROOFIT_TESTSTATISTICS_likelihood_builders diff --git a/roofit/roofitcore/inc/TestStatistics/optional_parameter_types.h b/roofit/roofitcore/inc/TestStatistics/optional_parameter_types.h new file mode 100644 index 0000000000000..8e48621327180 --- /dev/null +++ b/roofit/roofitcore/inc/TestStatistics/optional_parameter_types.h @@ -0,0 +1,51 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#ifndef ROOT_ROOFIT_TESTSTATISTICS_optional_parameter_types +#define ROOT_ROOFIT_TESTSTATISTICS_optional_parameter_types + +#include + +namespace RooFit { +namespace TestStatistics { + +// strongly named container types for use as optional parameters in the test statistics constructors + +/// Optional parameter used in buildLikelihood(), see documentation there. +struct ConstrainedParameters { + ConstrainedParameters(); + explicit ConstrainedParameters(const RooArgSet ¶meters); + RooArgSet set; +}; + +/// Optional parameter used in buildLikelihood(), see documentation there. +struct ExternalConstraints { + ExternalConstraints(); + explicit ExternalConstraints(const RooArgSet &constraints); + RooArgSet set; +}; + +/// Optional parameter used in buildLikelihood(), see documentation there. +struct GlobalObservables { + GlobalObservables(); + explicit GlobalObservables(const RooArgSet &observables); + RooArgSet set; +}; + +} +} + +#endif // ROOT_ROOFIT_TESTSTATISTICS_optional_parameter_types diff --git a/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx b/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx index 233f3759c0daf..ee578f60ce577 100644 --- a/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooAbsMinimizerFcn.cxx @@ -59,7 +59,7 @@ RooAbsMinimizerFcn::RooAbsMinimizerFcn(RooArgList paramList, RooMinimizer *conte for (unsigned int i = 0; i < _floatParamList->size(); ) { // Note: Counting loop, since removing from collection! const RooAbsArg* arg = (*_floatParamList).at(i); if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) { - oocoutW(_context,Minimization) << "RooMinimizerFcn::RooMinimizerFcn: removing parameter " + oocoutW(_context,Minimization) << "RooAbsMinimizerFcn::RooAbsMinimizerFcn: removing parameter " << arg->GetName() << " from list because it is not of type RooRealVar" << endl; _floatParamList->remove(*arg); } else { @@ -424,11 +424,11 @@ void RooAbsMinimizerFcn::printEvalErrors() const { std::ostringstream msg; if (_doEvalErrorWall) { - msg << "RooMinimizerFcn: Minimized function has error status." << endl + msg << "RooAbsMinimizerFcn: Minimized function has error status." << endl << "Returning maximum FCN so far (" << _maxFCN << ") to force MIGRAD to back out of this region. Error log follows.\n"; } else { - msg << "RooMinimizerFcn: Minimized function has error status but is ignored.\n"; + msg << "RooAbsMinimizerFcn: Minimized function has error status but is ignored.\n"; } msg << "Parameter values: " ; diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 7bcdab4c708c7..4b0706948728e 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -38,21 +38,9 @@ Various methods are available to control verbosity, profiling, automatic PDF optimization. **/ +#include "RooMinimizer.h" #include "RooFit.h" - -#include "TClass.h" - -#include -#include - -#include "TH2.h" -#include "TMarker.h" -#include "TGraph.h" -#include "Fit/FitConfig.h" -#include "TStopwatch.h" -#include "TMatrixDSym.h" - #include "RooArgSet.h" #include "RooArgList.h" #include "RooAbsReal.h" @@ -62,11 +50,21 @@ automatic PDF optimization. #include "RooSentinel.h" #include "RooMsgService.h" #include "RooPlot.h" - +#include "RooMinimizerFcn.h" +#include "RooGradMinimizerFcn.h" #include "RooFitResult.h" -#include "RooMinimizer.h" +#include "TClass.h" #include "Math/Minimizer.h" +#include "TH2.h" +#include "TMarker.h" +#include "TGraph.h" +#include "Fit/FitConfig.h" +#include "TStopwatch.h" +#include "TMatrixDSym.h" + +#include +#include #if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3) char* operator+( streampos&, char* ); @@ -251,9 +249,6 @@ void RooMinimizer::setOffsetting(Bool_t flag) //////////////////////////////////////////////////////////////////////////////// /// Choose the minimizer algorithm. -// forward declaration (necessary for avoiding circular dependency problems) -class RooGradMinimizerFcn; - void RooMinimizer::setMinimizerType(const char* type) { if (_fcnMode != FcnMode::classic && strcmp(type, "Minuit2") != 0) { @@ -329,6 +324,10 @@ bool RooMinimizer::fitFcn() const { ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); break; } + case FcnMode::generic_wrapper: { + ret = _theFitter->FitFCN(*dynamic_cast(_fcn)); + break; + } default: { throw std::logic_error("In RooMinimizer::fitFcn: _fcnMode has an unsupported value!"); } @@ -872,6 +871,9 @@ ROOT::Math::IMultiGenFunction* RooMinimizer::getMultiGenFcn() const case FcnMode::gradient: { return static_cast(dynamic_cast(_fcn)); } + case FcnMode::generic_wrapper: { + return static_cast(dynamic_cast(_fcn)); + } default: { throw std::logic_error("In RooMinimizer::getMultiGenFcn: _fcnMode has an unsupported value!"); } @@ -890,6 +892,9 @@ const RooAbsMinimizerFcn *RooMinimizer::fitterFcn() const case FcnMode::gradient: { return static_cast(dynamic_cast(getFitterMultiGenFcn())); } + case FcnMode::generic_wrapper: { + return static_cast(dynamic_cast(getFitterMultiGenFcn())); + } default: { throw std::logic_error("In RooMinimizer::fitterFcn: _fcnMode has an unsupported value!"); } diff --git a/roofit/roofitcore/src/RooNLLVar.cxx b/roofit/roofitcore/src/RooNLLVar.cxx index 540647bb463c2..92d73bc54858f 100644 --- a/roofit/roofitcore/src/RooNLLVar.cxx +++ b/roofit/roofitcore/src/RooNLLVar.cxx @@ -459,6 +459,16 @@ Double_t RooNLLVar::evaluatePartition(std::size_t firstEvent, std::size_t lastEv /// \param[in] lastEvent First event not to be processed. /// \return Tuple with (Kahan sum of probabilities, carry of kahan sum, sum of weights) RooNLLVar::ComputeResult RooNLLVar::computeBatched(std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent) const +{ + auto pdfClone = static_cast(_funcClone); + return computeBatchedFunc(pdfClone, _dataClone, _evalData, _normSet, _weightSq, stepSize, firstEvent, lastEvent); +} + +// static function, also used from TestStatistics::RooUnbinnedL +RooNLLVar::ComputeResult RooNLLVar::computeBatchedFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, + std::unique_ptr &evalData, + RooArgSet *normSet, bool weightSq, std::size_t stepSize, + std::size_t firstEvent, std::size_t lastEvent) { const auto nEvents = lastEvent - firstEvent; @@ -466,37 +476,35 @@ RooNLLVar::ComputeResult RooNLLVar::computeBatched(std::size_t stepSize, std::si throw std::invalid_argument(std::string("Error in ") + __FILE__ + ": Step size for batch computations can only be 1."); } - auto pdfClone = static_cast(_funcClone); - // Create a RunContext that will own the memory where computation results are stored. // Holding on to this struct in between function calls will make sure that the memory // is only allocated once. - if (!_evalData) { - _evalData.reset(new RooBatchCompute::RunContext); + if (!evalData) { + evalData.reset(new RooBatchCompute::RunContext); } - _evalData->clear(); - _dataClone->getBatches(*_evalData, firstEvent, nEvents); + evalData->clear(); + dataClone->getBatches(*evalData, firstEvent, nEvents); - auto results = pdfClone->getLogProbabilities(*_evalData, _normSet); + auto results = pdfClone->getLogProbabilities(*evalData, normSet); #ifdef ROOFIT_CHECK_CACHED_VALUES for (std::size_t evtNo = firstEvent; evtNo < std::min(lastEvent, firstEvent + 10); ++evtNo) { - _dataClone->get(evtNo); - if (_dataClone->weight() == 0.) // 0-weight events are not cached, so cannot compare against them. + dataClone->get(evtNo); + if (dataClone->weight() == 0.) // 0-weight events are not cached, so cannot compare against them. continue; - assert(_dataClone->valid()); + assert(dataClone->valid()); try { // Cross check results with strict tolerance and complain - BatchInterfaceAccessor::checkBatchComputation(*pdfClone, *_evalData, evtNo-firstEvent, _normSet, 1.E-13); + BatchInterfaceAccessor::checkBatchComputation(*pdfClone, *evalData, evtNo-firstEvent, normSet, 1.E-13); } catch (std::exception& e) { std::cerr << __FILE__ << ":" << __LINE__ << " ERROR when checking batch computation for event " << evtNo << ":\n" << e.what() << std::endl; // It becomes a real problem if it's very wrong. We fail in this case: try { - BatchInterfaceAccessor::checkBatchComputation(*pdfClone, *_evalData, evtNo-firstEvent, _normSet, 1.E-9); + BatchInterfaceAccessor::checkBatchComputation(*pdfClone, *evalData, evtNo-firstEvent, normSet, 1.E-9); } catch (std::exception& e2) { assert(false); } @@ -507,9 +515,9 @@ RooNLLVar::ComputeResult RooNLLVar::computeBatched(std::size_t stepSize, std::si // Compute sum of event weights. First check if we need squared weights - const RooSpan eventWeights = _dataClone->getWeightBatch(firstEvent, nEvents); + const RooSpan eventWeights = dataClone->getWeightBatch(firstEvent, nEvents); //Capture member for lambda: - const bool retrieveSquaredWeights = _weightSq; + const bool retrieveSquaredWeights = weightSq; auto retrieveWeight = [&eventWeights, retrieveSquaredWeights](std::size_t i) { return retrieveSquaredWeights ? eventWeights[i] * eventWeights[i] : eventWeights[i]; }; @@ -519,7 +527,7 @@ RooNLLVar::ComputeResult RooNLLVar::computeBatched(std::size_t stepSize, std::si double uniformSingleEventWeight{0.0}; double sumOfWeights; if (eventWeights.empty()) { - uniformSingleEventWeight = retrieveSquaredWeights ? _dataClone->weightSquared() : _dataClone->weight(); + uniformSingleEventWeight = retrieveSquaredWeights ? dataClone->weightSquared() : dataClone->weight(); sumOfWeights = nEvents * uniformSingleEventWeight; for (std::size_t i = 0; i < results.size(); ++i) { //CHECK_VECTORISE kahanProb.AddIndexed(-uniformSingleEventWeight * results[i], i); @@ -567,21 +575,28 @@ RooNLLVar::ComputeResult RooNLLVar::computeBatched(std::size_t stepSize, std::si RooNLLVar::ComputeResult RooNLLVar::computeScalar(std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent) const { auto pdfClone = static_cast(_funcClone); + return computeScalarFunc(pdfClone, _dataClone, _normSet, _weightSq, stepSize, firstEvent, lastEvent); +} +// static function, also used from TestStatistics::RooUnbinnedL +RooNLLVar::ComputeResult RooNLLVar::computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, + RooArgSet *normSet, bool weightSq, std::size_t stepSize, + std::size_t firstEvent, std::size_t lastEvent) +{ ROOT::Math::KahanSum kahanWeight; ROOT::Math::KahanSum kahanProb; RooNaNPacker packedNaN(0.f); for (auto i=firstEvent; iget(i) ; + dataClone->get(i) ; - if (!_dataClone->valid()) continue; + if (!dataClone->valid()) continue; - Double_t eventWeight = _dataClone->weight(); //FIXME + Double_t eventWeight = dataClone->weight(); //FIXME if (0. == eventWeight * eventWeight) continue ; - if (_weightSq) eventWeight = _dataClone->weightSquared() ; + if (weightSq) eventWeight = dataClone->weightSquared() ; - const double term = -eventWeight * pdfClone->getLogVal(_normSet); + const double term = -eventWeight * pdfClone->getLogVal(normSet); kahanWeight.Add(eventWeight); kahanProb.Add(term); diff --git a/roofit/roofitcore/src/TestStatistics/ConstantTermsOptimizer.cxx b/roofit/roofitcore/src/TestStatistics/ConstantTermsOptimizer.cxx new file mode 100644 index 0000000000000..ad35a9c900a55 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/ConstantTermsOptimizer.cxx @@ -0,0 +1,170 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include + +#include +#include // complete type for dynamic cast +#include +#include +#include + +namespace RooFit { +namespace TestStatistics { + +/** \class ConstantTermsOptimizer + * + * \brief Analyzes a function given a dataset/observables for constant terms and caches those in the dataset + * + * This optimizer should be used on a consistent combination of function (usually a pdf) and a dataset with observables. + * It then analyzes the function to find parts that can be precalculated because they are constant given the set of + * observables. These are cached inside the dataset and used in subsequent evaluations of the function on that dataset. + * The typical use case for this is inside likelihood minimization where many calls of the same pdf/dataset combination + * are made. \p norm_set must provide the normalization set of the function, which would typically be the set of + * observables in the dataset; this is used to make sure all object caches are created before analysis by evaluating the + * function on this set at the beginning of enableConstantTermsOptimization. + */ + +RooArgSet ConstantTermsOptimizer::requiredExtraObservables() +{ + // TODO: the RooAbsOptTestStatistics::requiredExtraObservables() call this code was copied + // from was overloaded for RooXYChi2Var only; implement different options when necessary + return RooArgSet(); +} + +void ConstantTermsOptimizer::enableConstantTermsOptimization(RooAbsReal *function, RooArgSet *norm_set, RooAbsData *dataset, + bool applyTrackingOpt) +{ + // Trigger create of all object caches now in nodes that have deferred object creation + // so that cache contents can be processed immediately + function->getVal(norm_set); + + // Apply tracking optimization here. Default strategy is to track components + // of RooAddPdfs and RooRealSumPdfs. If these components are a RooProdPdf + // or a RooProduct respectively, track the components of these products instead + // of the product term + RooArgSet trackNodes; + + // Add safety check here - applyTrackingOpt will only be applied if present + // dataset is constructed in terms of a RooVectorDataStore + if (applyTrackingOpt) { + if (!dynamic_cast(dataset->store())) { + oocoutW((TObject *)nullptr, Optimization) + << "enableConstantTermsOptimization(function: " << function->GetName() + << ", dataset: " << dataset->GetName() + << ") WARNING Cache-and-track optimization (Optimize level 2) is only available for datasets" + << " implemented in terms of RooVectorDataStore - ignoring this option for current dataset" << std::endl; + applyTrackingOpt = kFALSE; + } + } + + if (applyTrackingOpt) { + RooArgSet branches; + function->branchNodeServerList(&branches); + for (const auto arg : branches) { + arg->setCacheAndTrackHints(trackNodes); + } + // Do not set CacheAndTrack on constant expressions + auto constNodes = (RooArgSet *)trackNodes.selectByAttrib("Constant", kTRUE); + trackNodes.remove(*constNodes); + delete constNodes; + + // Set CacheAndTrack flag on all remaining nodes + trackNodes.setAttribAll("CacheAndTrack", kTRUE); + } + + // Find all nodes that depend exclusively on constant parameters + RooArgSet cached_nodes; + + function->findConstantNodes(*dataset->get(), cached_nodes); + + // Cache constant nodes with dataset - also cache entries corresponding to zero-weights in data when using + // BinnedLikelihood + // NOTE: we pass nullptr as cache-owner here, because we don't use the cacheOwner() anywhere in TestStatistics + // TODO: make sure this (nullptr) is always correct + dataset->cacheArgs(nullptr, cached_nodes, norm_set, !function->getAttribute("BinnedLikelihood")); + + // Put all cached nodes in AClean value caching mode so that their evaluate() is never called + for (const auto cacheArg : cached_nodes) { + cacheArg->setOperMode(RooAbsArg::AClean); + } + + std::unique_ptr constNodes {(RooArgSet *)cached_nodes.selectByAttrib("ConstantExpressionCached", kTRUE)}; + RooArgSet actualTrackNodes(cached_nodes); + actualTrackNodes.remove(*constNodes); + if (constNodes->getSize() > 0) { + if (constNodes->getSize() < 20) { + oocoutI((TObject*)nullptr, Minimization) + << " The following expressions have been identified as constant and will be precalculated and cached: " + << *constNodes << std::endl; + } else { + oocoutI((TObject*)nullptr, Minimization) << " A total of " << constNodes->getSize() + << " expressions have been identified as constant and will be precalculated and cached." + << std::endl; + } + } + if (actualTrackNodes.getSize() > 0) { + if (actualTrackNodes.getSize() < 20) { + oocoutI((TObject*)nullptr, Minimization) << " The following expressions will be evaluated in cache-and-track mode: " + << actualTrackNodes << std::endl; + } else { + oocoutI((TObject*)nullptr, Minimization) << " A total of " << constNodes->getSize() + << " expressions will be evaluated in cache-and-track-mode." << std::endl; + } + } + + // Disable reading of observables that are no longer used + dataset->optimizeReadingWithCaching(*function, cached_nodes, requiredExtraObservables()); +} + +void ConstantTermsOptimizer::disableConstantTermsOptimization(RooAbsReal *function, RooArgSet *norm_set, RooArgSet *observables, + RooAbsData *dataset) +{ + // Delete the cache + dataset->resetCache(); + + // Reactivate all tree branches + dataset->setArgStatus(*dataset->get(), kTRUE); + + // Reset all nodes to ADirty + optimizeCaching(function, norm_set, observables, dataset); + + // Disable propagation of dirty state flags for observables + dataset->setDirtyProp(kFALSE); + + // _cachedNodes.removeAll(); + + // _optimized = kFALSE; +} + +void ConstantTermsOptimizer::optimizeCaching(RooAbsReal *function, RooArgSet *norm_set, RooArgSet *observables, RooAbsData *dataset) +{ + // Trigger create of all object caches now in nodes that have deferred object creation + // so that cache contents can be processed immediately + function->getVal(norm_set); + + // Set value caching mode for all nodes that depend on any of the observables to ADirty + function->optimizeCacheMode(*observables); + + // Disable propagation of dirty state flags for observables + dataset->setDirtyProp(kFALSE); + + // Disable reading of observables that are not used + dataset->optimizeReadingWithCaching(*function, RooArgSet(), requiredExtraObservables()) ; +} + +} // namespace TestStatistics +} // namespace RooFit \ No newline at end of file diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx new file mode 100644 index 0000000000000..ae0aee0f42f69 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientWrapper.cxx @@ -0,0 +1,61 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include "TestStatistics/LikelihoodGradientWrapper.h" +#include "RooMinimizer.h" + +namespace RooFit { +namespace TestStatistics { + +/** \class LikelihoodGradientWrapper + * \brief Virtual base class for implementation of likelihood gradient calculation strategies + * + * This class provides the interface necessary for RooMinimizer (through MinuitFcnGrad) to get the likelihood gradient + * values it needs for fitting the pdf to the data. The strategy by which these values are obtained is up to the + * implementer of this class. Its intended purpose was mainly to allow for parallel calculation strategies. + * + * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived + * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). + */ + +/* + * \param[in] likelihood Shared pointer to the likelihood that must be evaluated + * \param[in] calculation_is_clean Shared pointer to the object that keeps track of what has been evaluated for the current parameter set provided by Minuit. This information can be used by different calculators, so must be shared between them. + * \param[in] minimizer Raw pointer to the minimizer that owns the MinuitFcnGrad object that owns this wrapper object. + */ +LikelihoodGradientWrapper::LikelihoodGradientWrapper(std::shared_ptr likelihood, + std::shared_ptr calculation_is_clean, + std::size_t /*N_dim*/, RooMinimizer *minimizer) + : likelihood_(std::move(likelihood)), minimizer_(minimizer), calculation_is_clean_(std::move(calculation_is_clean)) +{ + // Note to future maintainers: take care when storing the minimizer_fcn pointer. The + // RooAbsMinimizerFcn subclasses may get cloned inside MINUIT, which means the pointer + // should also somehow be updated in this class. +} + +void LikelihoodGradientWrapper::synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions & /*options*/) {} + +void LikelihoodGradientWrapper::synchronizeParameterSettings( + const std::vector ¶meter_settings) +{ + synchronizeParameterSettings(minimizer_->getMultiGenFcn(), parameter_settings); +} + +void LikelihoodGradientWrapper::updateMinuitInternalParameterValues(const std::vector& /*minuit_internal_x*/) {} +void LikelihoodGradientWrapper::updateMinuitExternalParameterValues(const std::vector& /*minuit_external_x*/) {} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodSerial.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodSerial.cxx new file mode 100644 index 0000000000000..971338139958c --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodSerial.cxx @@ -0,0 +1,104 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + + +#include +#include +#include +#include +#include +#include +#include "RooRealVar.h" + +namespace RooFit { +namespace TestStatistics { + +/** \class LikelihoodSerial + * \brief Serial likelihood calculation strategy implementation + * + * This class serves as a baseline reference implementation of the LikelihoodWrapper. It reimplements the previous + * RooNLLVar "BulkPartition" single CPU strategy in the new RooFit::TestStatistics framework. + * + * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived + * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). + */ + +LikelihoodSerial::LikelihoodSerial(std::shared_ptr likelihood, std::shared_ptr calculation_is_clean) + : LikelihoodWrapper(std::move(likelihood), std::move(calculation_is_clean)) +{ + initVars(); + // determine likelihood type + if (dynamic_cast(likelihood_.get()) != nullptr) { + likelihood_type = LikelihoodType::unbinned; + } else if (dynamic_cast(likelihood_.get()) != nullptr) { + likelihood_type = LikelihoodType::binned; + } else if (dynamic_cast(likelihood_.get()) != nullptr) { + likelihood_type = LikelihoodType::sum; + } else if (dynamic_cast(likelihood_.get()) != nullptr) { + likelihood_type = LikelihoodType::subsidiary; + } else { + throw std::logic_error("in LikelihoodSerial constructor: _likelihood is not of a valid subclass!"); + } + // Note to future maintainers: take care when storing the minimizer_fcn pointer. The + // RooAbsMinimizerFcn subclasses may get cloned inside MINUIT, which means the pointer + // should also somehow be updated in this class. +} + +/// \brief Helper function for the constuctor. +/// +/// This is a separate function (instead of just in ctor) for historical reasons. +/// Its predecessor RooRealMPFE::initVars() was used from multiple ctors, but also +/// from RooRealMPFE::constOptimizeTestStatistic at the end, which makes sense, +/// because it might change the set of variables. We may at some point want to do +/// this here as well. +void LikelihoodSerial::initVars() +{ + // Empty current lists + _vars.removeAll(); + _saveVars.removeAll(); + + // Retrieve non-constant parameters + auto vars = std::make_unique(*likelihood_->getParameters()); + + RooArgList varList(*vars); + + // Save in lists + _vars.add(varList); + _saveVars.addClone(varList); +} + +void LikelihoodSerial::evaluate() { + switch (likelihood_type) { + case LikelihoodType::unbinned: + case LikelihoodType::binned: { + result = likelihood_->evaluatePartition({0, 1}, 0, 0); + break; + } + case LikelihoodType::sum: { + result = likelihood_->evaluatePartition({0, 1}, 0, likelihood_->getNComponents()); + break; + } + default: { + throw std::logic_error("in LikelihoodSerial::evaluate_task: likelihood types other than binned, unbinned and simultaneous not yet implemented!"); + break; + } + } + + result = applyOffsetting(result); +} + +} // namespace TestStatistics +} // namespace RooFit \ No newline at end of file diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodWrapper.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodWrapper.cxx new file mode 100644 index 0000000000000..4370839099ad4 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodWrapper.cxx @@ -0,0 +1,151 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include + +#include // need complete type for likelihood->... +#include +#include +#include // need complete type for dynamic cast + +#include "RooMsgService.h" + +namespace RooFit { +namespace TestStatistics { + +/** \class LikelihoodWrapper + * \brief Virtual base class for implementation of likelihood calculation strategies + * + * This class provides the interface necessary for RooMinimizer (through MinuitFcnGrad) to get the likelihood values it + * needs for fitting the pdf to the data. The strategy by which these values are obtained is up to the implementer of + * this class. Its intended purpose was mainly to allow for parallel calculation strategies, but serial strategies are + * possible too, as illustrated in LikelihoodSerial. + * + * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived + * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). + */ + +/* + * \param[in] likelihood Shared pointer to the likelihood that must be evaluated + * \param[in] calculation_is_clean Shared pointer to the object that keeps track of what has been evaluated for the current parameter set provided by Minuit. This information can be used by different calculators, so must be shared between them. + */ +LikelihoodWrapper::LikelihoodWrapper(std::shared_ptr likelihood, + std::shared_ptr calculation_is_clean) + : likelihood_(std::move(likelihood)), + calculation_is_clean_(std::move(calculation_is_clean)) +{ + // Note to future maintainers: take care when storing the minimizer_fcn pointer. The + // RooAbsMinimizerFcn subclasses may get cloned inside MINUIT, which means the pointer + // should also somehow be updated in this class. +} + +void LikelihoodWrapper::synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions & /*options*/) {} + +void LikelihoodWrapper::synchronizeParameterSettings( + const std::vector & /*parameter_settings*/) +{ +} + +void LikelihoodWrapper::constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt) +{ + likelihood_->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); +} + +double LikelihoodWrapper::defaultErrorLevel() const +{ + return likelihood_->defaultErrorLevel(); +} +std::string LikelihoodWrapper::GetName() const +{ + return likelihood_->GetName(); +} +std::string LikelihoodWrapper::GetTitle() const +{ + return likelihood_->GetTitle(); +} + +void LikelihoodWrapper::enableOffsetting(bool flag) +{ + do_offset_ = flag; + // Clear offset if feature is disabled so that it is recalculated next time it is enabled + if (!do_offset_) { + offset_ = {}; + } +} + +void LikelihoodWrapper::setOffsettingMode(OffsettingMode mode) +{ + offsetting_mode_ = mode; + if (isOffsetting()) { + oocoutI(static_cast(nullptr), Minimization) << "LikelihoodWrapper::setOffsettingMode(" << GetName() << "): changed offsetting mode while offsetting was enabled; resetting offset values" << std::endl; + offset_ = {}; + } +} + +ROOT::Math::KahanSum LikelihoodWrapper::applyOffsetting(ROOT::Math::KahanSum current_value) +{ + if (do_offset_) { + + // If no offset is stored enable this feature now + if (offset_ == 0 && current_value != 0) { + offset_ = current_value; + if (offsetting_mode_ == OffsettingMode::legacy) { + auto sum_likelihood = dynamic_cast(likelihood_.get()); + if (sum_likelihood != nullptr) { + auto subsidiary_value = sum_likelihood->getSubsidiaryValue(); + // "undo" the addition of the subsidiary value to emulate legacy behavior + offset_ -= subsidiary_value; + // manually calculate result with zero carry, again to emulate legacy behavior + return {current_value.Result() - offset_.Result()}; + } + } + oocoutI(static_cast(nullptr), Minimization) + << "LikelihoodWrapper::applyOffsetting(" << GetName() << "): Likelihood offset now set to " << offset_ + << std::endl; + } + + return current_value - offset_; + } else { + return current_value; + } +} + +/// When calculating an unbinned likelihood with square weights applied, a different offset +/// is necessary. Similar situations may ask for a separate offset as well. This function +/// switches between the two sets of offset values. +void LikelihoodWrapper::swapOffsets() +{ + std::swap(offset_, offset_save_); +} + +void LikelihoodWrapper::setApplyWeightSquared(bool flag) +{ + RooUnbinnedL *unbinned_likelihood = dynamic_cast(likelihood_.get()); + if (unbinned_likelihood == nullptr) { + throw std::logic_error("LikelihoodWrapper::setApplyWeightSquared can only be used on unbinned likelihoods, but the wrapped likelihood_ member is not a RooUnbinnedL!"); + } + bool flag_was_changed = unbinned_likelihood->setApplyWeightSquared(flag); + + if (flag_was_changed) { + swapOffsets(); + } +} + +void LikelihoodWrapper::updateMinuitInternalParameterValues(const std::vector& /*minuit_internal_x*/) {} +void LikelihoodWrapper::updateMinuitExternalParameterValues(const std::vector& /*minuit_external_x*/) {} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx new file mode 100644 index 0000000000000..46a9e68a3b2b5 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -0,0 +1,220 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include "TestStatistics/MinuitFcnGrad.h" + +#include "RooMinimizer.h" +#include "RooMsgService.h" +#include "RooAbsPdf.h" + +#include // std::setprecision + +namespace RooFit { +namespace TestStatistics { + +/** \class MinuitFcnGrad + * + * \brief Minuit-RooMinimizer interface which synchronizes parameter data and coordinates evaluation of likelihood (gradient) values + * + * This class provides an interface between RooFit and Minuit. It synchronizes parameter values from Minuit, calls + * calculator classes to evaluate likelihood and likelihood gradient values and returns them to Minuit. The Wrapper + * objects do the actual calculations. These are constructed inside the MinuitFcnGrad constructor using the RooAbsL + * likelihood passed in to the constructor, usually directly from RooMinimizer, with which this class is intimately + * coupled, being a RooAbsMinimizerFcn implementation. MinuitFcnGrad inherits from ROOT::Math::IMultiGradFunction as + * well, which allows it to be used as the FCN and GRAD parameters Minuit expects. + * + * \note The class is not intended for use by end-users. We recommend to either use RooMinimizer with a RooAbsL derived + * likelihood object, or to use a higher level entry point like RooAbsPdf::fitTo() or RooAbsPdf::createNLL(). + */ + +double MinuitFcnGrad::DoEval(const double *x) const +{ + Bool_t parameters_changed = syncParameterValuesFromMinuitCalls(x, false); + + // Calculate the function for these parameters +// RooAbsReal::setHideOffset(kFALSE); + likelihood->evaluate(); + double fvalue = likelihood->getResult().Sum(); + calculation_is_clean->likelihood = true; +// RooAbsReal::setHideOffset(kTRUE); + + if (!parameters_changed) { + return fvalue; + } + + if (!std::isfinite(fvalue) || RooAbsReal::numEvalErrors() > 0 || fvalue > 1e30) { + + if (_printEvalErrors >= 0) { + + if (_doEvalErrorWall) { + oocoutW(static_cast(nullptr), Eval) + << "RooGradMinimizerFcn: Minimized function has error status." << std::endl + << "Returning maximum FCN so far (" << _maxFCN + << ") to force MIGRAD to back out of this region. Error log follows" << std::endl; + } else { + oocoutW(static_cast(nullptr), Eval) + << "RooGradMinimizerFcn: Minimized function has error status but is ignored" << std::endl; + } + + Bool_t first(kTRUE); + ooccoutW(static_cast(nullptr), Eval) << "Parameter values: "; + for (const auto rooAbsArg : *_floatParamList) { + auto var = static_cast(rooAbsArg); + if (first) { + first = kFALSE; + } else { + ooccoutW(static_cast(nullptr), Eval) << ", "; + } + ooccoutW(static_cast(nullptr), Eval) << var->GetName() << "=" << var->getVal(); + } + ooccoutW(static_cast(nullptr), Eval) << std::endl; + + RooAbsReal::printEvalErrors(ooccoutW(static_cast(nullptr), Eval), _printEvalErrors); + ooccoutW(static_cast(nullptr), Eval) << std::endl; + } + + if (_doEvalErrorWall) { + fvalue = _maxFCN + 1; + } + + RooAbsReal::clearEvalErrorLog(); + _numBadNLL++; + } else if (fvalue > _maxFCN) { + _maxFCN = fvalue; + } + + // Optional logging + if (_verbose) { + std::cout << "\nprevFCN" << (likelihood->isOffsetting() ? "-offset" : "") << " = " << std::setprecision(10) + << fvalue << std::setprecision(4) << " "; + std::cout.flush(); + } + + _evalCounter++; + return fvalue; +} + +/// Minuit calls (via FcnAdapters etc) DoEval or Gradient with a set of parameters x. +/// This function syncs these values to the proper places in RooFit. +/// +/// The first twist, and reason this function is more complicated than one may imagine, is that Minuit internally uses a +/// transformed parameter space to account for parameter boundaries. Whether we receive these Minuit "internal" +/// parameter values or "regular"/untransformed RooFit parameter space values depends on the situation. +/// - The values that arrive here via DoEval are always "normal" parameter values, since Minuit transforms these +/// back into regular space before passing to DoEval (see MnUserFcn::operator() which wraps the Fcn(Gradient)Base +/// in ModularFunctionMinimizer::Minimize and is used for direct function calls from that point on in the minimizer). +/// These can thus always be safely synced with this function's RooFit parameters using SetPdfParamVal. +/// - The values that arrive here via Gradient will be in internal coordinates if that is +/// what this class expects, and indeed this is the case for MinuitFcnGrad's current implementation. This is +/// communicated to Minuit via MinuitFcnGrad::returnsInMinuit2ParameterSpace. Inside Minuit, that function determines +/// whether this class's gradient calculator is wrapped inside a AnalyticalGradientCalculator, to which Minuit passes +/// "external" parameter values, or as an ExternalInternalGradientCalculator, which gets "internal" parameter values. +/// Long story short: when MinuitFcnGrad::returnsInMinuit2ParameterSpace() returns true, Minuit will pass "internal" +/// values to Gradient. These cannot be synced with this function's RooFit parameters using +/// SetPdfParamVal, unless a manual transformation step is performed in advance. However, they do need to be passed +/// on to the gradient calculator, since indeed we expect values there to be in "internal" space. However, this is +/// calculator dependent. Note that in the current MinuitFcnGrad implementation we do not actually allow for +/// calculators in "external" (i.e. regular RooFit parameter space) values, since +/// MinuitFcnGrad::returnsInMinuit2ParameterSpace is hardcoded to true. This should in a future version be changed so +/// that the calculator (the wrapper) is queried for this information. +/// Because some gradient calculators may also use the regular RooFit parameters (e.g. for calculating the likelihood's +/// value itself), this information is also passed on to the gradient wrapper. Vice versa, when updated "internal" +/// parameters are passed to Gradient, the likelihood may be affected as well. Even though a +/// transformation from internal to "external" may be necessary before the values can be used, the likelihood can at +/// least log that its parameter values are possibly no longer in sync with those of the gradient. +/// +/// The second twist is that the Minuit external parameters may still be different from the ones used in RooFit. This +/// happens when Minuit tries out values that lay outside the RooFit parameter's range(s). RooFit's setVal (called +/// inside SetPdfParamVal) then clips the RooAbsArg's value to one of the range limits, instead of setting it to the +/// value Minuit intended. When this happens, i.e. syncParameterValuesFromMinuitCalls is called with +/// minuit_internal = false and the values do not match the previous values stored in minuit_internal_x_ *but* the +/// values after SetPdfParamVal did not get set to the intended value, the minuit_internal_roofit_x_mismatch_ flag is +/// set. This information can be used by calculators, if desired, for instance when a calculator does not want to make +/// use of the range information in the RooAbsArg parameters. +bool MinuitFcnGrad::syncParameterValuesFromMinuitCalls(const double *x, bool minuit_internal) const +{ + bool a_parameter_has_been_updated = false; + if (minuit_internal) { + if (!returnsInMinuit2ParameterSpace()) { + throw std::logic_error("Updating Minuit-internal parameters only makes sense for (gradient) calculators that are defined in Minuit-internal parameter space."); + } + + for (std::size_t ix = 0; ix < NDim(); ++ix) { + bool parameter_changed = (x[ix] != minuit_internal_x_[ix]); + if (parameter_changed) { + minuit_internal_x_[ix] = x[ix]; + } + a_parameter_has_been_updated |= parameter_changed; + } + + if(a_parameter_has_been_updated) { + calculation_is_clean->set_all(false); + likelihood->updateMinuitInternalParameterValues(minuit_internal_x_); + gradient->updateMinuitInternalParameterValues(minuit_internal_x_); + } + } else { + bool a_parameter_is_mismatched = false; + + for (std::size_t ix = 0; ix < NDim(); ++ix) { + // Note: the return value of SetPdfParamVal does not always mean that the parameter's value in the RooAbsReal changed since last + // time! If the value was out of range bin, setVal was still called, but the value was not updated. + SetPdfParamVal(ix, x[ix]); + minuit_external_x_[ix] = x[ix]; + // The above is why we need minuit_external_x_. The minuit_external_x_ vector can also be passed to + // LikelihoodWrappers, if needed, but typically they will make use of the RooFit parameters directly. However, + // we log in the flag below whether they are different so that calculators can use this information. + bool parameter_changed = (x[ix] != minuit_external_x_[ix]); + a_parameter_has_been_updated |= parameter_changed; + a_parameter_is_mismatched |= (((RooRealVar *)_floatParamList->at(ix))->getVal() != minuit_external_x_[ix]); + } + + minuit_internal_roofit_x_mismatch_ = a_parameter_is_mismatched; + + if(a_parameter_has_been_updated) { + calculation_is_clean->set_all(false); + likelihood->updateMinuitExternalParameterValues(minuit_external_x_); + gradient->updateMinuitExternalParameterValues(minuit_external_x_); + } + } + return a_parameter_has_been_updated; +} + + +void MinuitFcnGrad::Gradient(const double *x, double *grad) const +{ + syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); + gradient->fillGradient(grad); +} + +double MinuitFcnGrad::DoDerivative(const double * /*x*/, unsigned int /*icoord*/) const +{ + throw std::runtime_error("MinuitFcnGrad::DoDerivative is not implemented, please use Gradient instead."); +} + +Bool_t +MinuitFcnGrad::Synchronize(std::vector ¶meters, Bool_t optConst, Bool_t verbose) +{ + Bool_t returnee = synchronizeParameterSettings(parameters, optConst, verbose); + likelihood->synchronizeParameterSettings(parameters); + gradient->synchronizeParameterSettings(parameters); + + likelihood->synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); + gradient->synchronizeWithMinimizer(_context->fitter()->Config().MinimizerOptions()); + return returnee; +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx b/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx new file mode 100644 index 0000000000000..831335013b675 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooAbsL.cxx @@ -0,0 +1,254 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include +#include +#include "RooAbsPdf.h" +#include "RooAbsData.h" + +// 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 pdf, std::shared_ptr 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((TObject *)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(std::shared_ptr(static_cast(in.pdf->cloneTree())), + std::shared_ptr(static_cast(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(nullptr), inpdf}, {std::shared_ptr(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(other.pdf_->cloneTree())); + data_.reset(static_cast(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 origParams{inpdf.getParameters(indata)}; + pdf_->recursiveRedirectServers(*origParams); + + // Store normalization set + normSet_.reset((RooArgSet *)indata.get()->snapshot(kFALSE)); + + // Expand list of observables with any observables used in parameterized ranges + for (const auto realDep : *_funcObsSet) { + auto realDepRLV = dynamic_cast(realDep); + if (realDepRLV && realDepRLV->isDerived()) { + RooArgSet tmp2; + realDepRLV->leafNodeServerList(&tmp2, 0, kTRUE); + _funcObsSet->add(tmp2, kTRUE); + } + } + + // ****************************************************************** + // *** 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(arg); + if (!realReal) { + continue; + } + auto datReal = dynamic_cast(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((TObject *)0, 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((TObject *)0, 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 + + // 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(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(kFALSE); + + // Disable reading of observables that are not used + data_->optimizeReadingWithCaching(*pdf_, RooArgSet(), RooArgSet()); +} + +RooArgSet *RooAbsL::getParameters() +{ + auto ding = pdf_->getParameters(*data_); + return ding; +} + +void RooAbsL::constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt) +{ + // to be further implemented, this is just a first test implementation + if (opcode == RooAbsArg::Activate) { + ConstantTermsOptimizer::enableConstantTermsOptimization(pdf_.get(), normSet_.get(), data_.get(), doAlsoTrackingOpt); + } +} + +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(data_->numEntries()); +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooBinnedL.cxx b/roofit/roofitcore/src/TestStatistics/RooBinnedL.cxx new file mode 100644 index 0000000000000..b8f18bdc23c29 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooBinnedL.cxx @@ -0,0 +1,151 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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 RooBinnedL.cxx +\class RooBinnedL +\ingroup Roofitcore + +Class RooBinnedL implements a -log(likelihood) calculation from a dataset +(assumed to be binned) and a PDF. The NLL is calculated as +\f + \sum_\mathrm{data} -\log( \mathrm{pdf}(x_\mathrm{data})) +\f] +In extended mode, a +\f$ N_\mathrm{expect} - N_\mathrm{observed}*log(N_\mathrm{expect}) \f$ term is added. +**/ + +#include + +#include "RooAbsData.h" +#include "RooAbsPdf.h" +#include "RooAbsDataStore.h" +#include "RooRealSumPdf.h" +#include "RooRealVar.h" + +#include "TMath.h" +#include "Math/Util.h" // KahanSum + +namespace RooFit { +namespace TestStatistics { + +RooBinnedL::RooBinnedL(RooAbsPdf *pdf, RooAbsData *data) + : RooAbsL(RooAbsL::ClonePdfData{pdf, data}, data->numEntries(), 1) +{ + // pdf must be a RooRealSumPdf representing a yield vector for a binned likelihood calculation + if (!dynamic_cast(pdf)) { + throw std::logic_error("RooBinnedL can only be created from pdf of type RooRealSumPdf!"); + } + + // Retrieve and cache bin widths needed to convert unnormalized binned pdf values back to yields + + // The Active label will disable pdf integral calculations + pdf->setAttribute("BinnedLikelihoodActive"); + + RooArgSet *obs = pdf->getObservables(data); + if (obs->getSize() != 1) { + throw std::logic_error( + "RooBinnedL can only be created from combination of pdf and data which has exactly one observable!"); + } else { + RooRealVar *var = (RooRealVar *)obs->first(); + std::list *boundaries = pdf->binBoundaries(*var, var->getMin(), var->getMax()); + std::list::iterator biter = boundaries->begin(); + _binw.resize(boundaries->size() - 1); + Double_t lastBound = (*biter); + ++biter; + int ibin = 0; + while (biter != boundaries->end()) { + _binw[ibin] = (*biter) - lastBound; + lastBound = (*biter); + ibin++; + ++biter; + } + } +} + +////////////////////////////////////////////////////////////////////////////////// +/// Calculate and return likelihood on subset of data from firstEvent to lastEvent +/// processed with a step size of 'stepSize'. If this an extended likelihood and +/// and the zero event is processed the extended term is added to the return +/// likelihood. +// +ROOT::Math::KahanSum +RooBinnedL::evaluatePartition(Section bins, std::size_t /*components_begin*/, std::size_t /*components_end*/) +{ + // Throughout the calculation, we use Kahan's algorithm for summing to + // prevent loss of precision - this is a factor four more expensive than + // straight addition, but since evaluating the PDF is usually much more + // expensive than that, we tolerate the additional cost... + ROOT::Math::KahanSum result; + +// data->store()->recalculateCache(_projDeps, firstEvent, lastEvent, stepSize, (_binnedPdf?kFALSE:kTRUE)); + // TODO: check when we might need _projDeps (it seems to be mostly empty); ties in with TODO below + data_->store()->recalculateCache(nullptr, bins.begin(N_events_), bins.end(N_events_), 1, kFALSE); + + ROOT::Math::KahanSum sumWeight; + + for (std::size_t i = bins.begin(N_events_); i < bins.end(N_events_); ++i) { + + data_->get(i); + + if (!data_->valid()) + continue; + + Double_t eventWeight = data_->weight(); + + // Calculate log(Poisson(N|mu) for this bin + Double_t N = eventWeight; + Double_t mu = pdf_->getVal() * _binw[i]; + + if (mu <= 0 && N > 0) { + + // Catch error condition: data present where zero events are predicted +// logEvalError(Form("Observed %f events in bin %d with zero event yield", N, i)); + // TODO: check if using regular stream vs logEvalError error gathering is ok + oocoutI(static_cast(nullptr), Minimization) + << "Observed " << N << " events in bin " << i << " with zero event yield" << std::endl; + + } else if (fabs(mu) < 1e-10 && fabs(N) < 1e-10) { + + // Special handling of this case since log(Poisson(0,0)=0 but can't be calculated with usual log-formula + // since log(mu)=0. No update of result is required since term=0. + + } else { + + Double_t term = -1 * (-mu + N * log(mu) - TMath::LnGamma(N + 1)); + + sumWeight += eventWeight; + result += term; + } + } + + // If part of simultaneous PDF normalize probability over + // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n) + if (sim_count_ > 1) { + result += sumWeight * log(1.0 * sim_count_); + } + + // At the end of the first full calculation, wire the caches + if (_first) { + _first = false; + pdf_->wireAllCaches(); + } + + return result; +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooRealL.cxx b/roofit/roofitcore/src/TestStatistics/RooRealL.cxx new file mode 100644 index 0000000000000..6afdfedd7293f --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooRealL.cxx @@ -0,0 +1,61 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include +#include + +namespace RooFit { +namespace TestStatistics { + +/** \class RooRealL + * \ingroup Roofitcore + * + * \brief RooAbsReal that wraps RooAbsL likelihoods for use in RooFit outside of the RooMinimizer context + * + * This class provides a simple wrapper to evaluate RooAbsL derived likelihood objects like a regular RooFit real value. + * Whereas the RooAbsL objects are meant to be used within the context of minimization, RooRealL can be used in any + * RooFit context, like plotting. The value can be accessed through getVal(), like with other RooFit real variables. + **/ + +RooRealL::RooRealL(const char *name, const char *title, std::shared_ptr likelihood) + : RooAbsReal(name, title), likelihood_(std::move(likelihood)), + vars_proxy_("varsProxy", "proxy set of parameters", this) +{ + vars_proxy_.add(*likelihood_->getParameters()); +} + +RooRealL::RooRealL(const RooRealL &other, const char *name) + : RooAbsReal(other, name), likelihood_(other.likelihood_), vars_proxy_("varsProxy", "proxy set of parameters", this) +{ + vars_proxy_.add(*likelihood_->getParameters()); +} + +Double_t RooRealL::evaluate() const +{ + // Evaluate as straight FUNC + std::size_t last_component = likelihood_->getNComponents(); + + auto ret_kahan = likelihood_->evaluatePartition({0, 1}, 0, last_component); + + const Double_t norm = globalNormalization(); + double ret = ret_kahan.Sum() / norm; + eval_carry = ret_kahan.Carry() / norm; + + return ret; +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooSubsidiaryL.cxx b/roofit/roofitcore/src/TestStatistics/RooSubsidiaryL.cxx new file mode 100644 index 0000000000000..a4e78e1f78731 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooSubsidiaryL.cxx @@ -0,0 +1,84 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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 RooSubsidiaryL.cxx +\class RooSubsidiaryL +\ingroup Roofitcore + +\brief RooSubsidiaryL calculates the sum of the -(log) likelihoods of a set of RooAbsPdf objects that represent +subsidiary or constraint functions. + +This class is used to gather all subsidiary PDF terms from the component PDFs of RooSumL likelihoods and calculate the +composite -log(L). Such subsidiary terms can be marked using RooFit::Constrain arguments to RooAbsPdf::fitTo() or +RooAbsPdf::createNLL(). + +Subsidiary terms are summed separately for increased numerical stability, since these terms are often small and cause +numerical variances in their original PDFs, whereas by summing as one separate subsidiary collective term, it is +numerically very stable. + +\note When a subsidiary PDF is part of multiple component PDFs, it will only be summed once in this class! This doesn't +change the derivative of the log likelihood (which is what matters in fitting the likelihood), but does change the +value of the (log-)likelihood itself. +**/ + +#include +#include // for dynamic cast +#include + +#include "Math/Util.h" // KahanSum + +namespace RooFit { +namespace TestStatistics { + +RooSubsidiaryL::RooSubsidiaryL(const std::string &parent_pdf_name, const RooArgSet &pdfs, + const RooArgSet ¶meter_set) + : RooAbsL(nullptr, nullptr, 0, 0, RooAbsL::Extended::No), parent_pdf_name_(parent_pdf_name) +{ + for (const auto comp : pdfs) { + if (!dynamic_cast(comp)) { + oocoutE((TObject *)0, InputArguments) << "RooSubsidiaryL::ctor(" << GetName() << ") ERROR: component " + << comp->GetName() << " is not of type RooAbsPdf" << std::endl; + RooErrorHandler::softAbort(); + } + subsidiary_pdfs_.add(*comp); + } + parameter_set_.add(parameter_set); +} + +ROOT::Math::KahanSum RooSubsidiaryL::evaluatePartition(RooAbsL::Section events, + std::size_t /*components_begin*/, + std::size_t /*components_end*/) +{ + if (events.begin_fraction != 0 || events.end_fraction != 1) { + oocoutW((TObject *)0, InputArguments) << "RooSubsidiaryL::evaluatePartition can only calculate everything, so " + "section should be {0,1}, but it's not!" + << std::endl; + } + + ROOT::Math::KahanSum sum; + + for (const auto comp : subsidiary_pdfs_) { + sum += -((RooAbsPdf *)comp)->getLogVal(¶meter_set_); + } + + return sum; +} + +void RooSubsidiaryL::constOptimizeTestStatistic(RooAbsArg::ConstOpCode /*opcode*/, bool /*doAlsoTrackingOpt*/) {} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooSumL.cxx b/roofit/roofitcore/src/TestStatistics/RooSumL.cxx new file mode 100644 index 0000000000000..698559aca9932 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooSumL.cxx @@ -0,0 +1,125 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include +#include +#include + +#include // min, max + +namespace RooFit { +namespace TestStatistics { + +/** \class RooSumL + * \ingroup Roofitcore + * + * \brief Likelihood class that sums over multiple -log components + * + * The likelihood is often a product of components, for instance when fitting simultaneous pdfs, but also when using + * subsidiary pdfs. Hence, the negative log likelihood that we, in fact, calculate is often a sum over these components. + * This sum is implemented by this class. + **/ + +/// \param[in] pdf Raw pointer to the pdf; will not be cloned in this object. +/// \param[in] data Raw pointer to the dataset; will not be cloned in this object. +/// \param[in] components The component likelihoods. +/// \param extended Set extended term calculation on, off or use Extended::Auto to determine automatically based on the +/// pdf whether to activate or not. \warning components must be passed with std::move, otherwise it cannot be moved into +/// the RooSumL because of the unique_ptr! \note The number of events in RooSumL is that of the full dataset. Components +/// will have their own number of events that may be more relevant. +RooSumL::RooSumL(RooAbsPdf *pdf, RooAbsData *data, std::vector> components, + RooAbsL::Extended extended) + : RooAbsL(pdf, data, data->numEntries(), components.size(), extended), components_(std::move(components)) +{ +} +// Developer note on the std::move() warning above: +// +// The point here was that you don't want to clone RooAbsL's too much, because they contain clones of the pdf and +// dataset that may have been mangled for optimization. You probably don't want to be doing that all the time, although +// it is a premature optimization, since we haven't timed its impact. That is the motivation behind using unique_ptrs +// for the components. The way the classes are built, the RooSumL doesn't care about what components it gets, so by +// definition it cannot create them internally, so they have to be passed in somehow. Forcing the user to call the +// function with a std::move is a way to make them fully realize that their local components will be destroyed and the +// contents moved into the RooSumL. +// +// We could change the type to an rvalue reference to make it clearer from the compiler error that std::move is +// necessary, instead of the obscure error that you get now. Compare the compiler error messages from these two types: +// +//#include +//#include +//#include +// +// struct Clear { +// Clear(std::vector>&& vec) : vec_(std::move(vec)) { +// printf("number is %d", *vec_[0]); +// } +// +// std::vector> vec_; +//}; +// +// struct Obscure { +// Obscure(std::vector> vec) : vec_(std::move(vec)) { +// printf("number is %d", *vec_[0]); +// } +// +// std::vector> vec_; +//}; +// +// int main() { +// std::vector> vec; +// vec.emplace_back(new int(4)); +// Clear thing(vec); +// Obscure thingy(vec); +//} + +/// \note Compared to the RooAbsTestStatistic implementation that this was taken from, we leave out Hybrid and +/// SimComponents interleaving support here. This should be implemented by a calculator (i.e. LikelihoodWrapper or +/// LikelihoodGradientWrapper derived class), if desired. +ROOT::Math::KahanSum +RooSumL::evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) +{ + // Evaluate specified range of owned GOF objects + ROOT::Math::KahanSum ret; + + // from RooAbsOptTestStatistic::combinedValue (which is virtual, so could be different for non-RooNLLVar!): + for (std::size_t ix = components_begin; ix < components_end; ++ix) { + ret += components_[ix]->evaluatePartition(events, 0, 0); + } + + return ret; +} + +/// \note This function assumes there is only one subsidiary component. +ROOT::Math::KahanSum RooSumL::getSubsidiaryValue() +{ + // iterate in reverse, because the subsidiary component is usually at the end: + for (auto component = components_.rbegin(); component != components_.rend(); ++component) { + if (dynamic_cast((*component).get()) != nullptr) { + return (*component)->evaluatePartition({0, 1}, 0, 0); + } + } + return {}; +} + +void RooSumL::constOptimizeTestStatistic(RooAbsArg::ConstOpCode opcode, bool doAlsoTrackingOpt) +{ + for (auto &component : components_) { + component->constOptimizeTestStatistic(opcode, doAlsoTrackingOpt); + } +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/RooUnbinnedL.cxx b/roofit/roofitcore/src/TestStatistics/RooUnbinnedL.cxx new file mode 100644 index 0000000000000..dbea172346032 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/RooUnbinnedL.cxx @@ -0,0 +1,171 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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 RooUnbinnedL.cxx +\class RooUnbinnedL +\ingroup Roofitcore + +Class RooUnbinnedL implements a -log(likelihood) calculation from a dataset +(assumed to be unbinned) and a PDF. The NLL is calculated as +\f[ + \sum_\mathrm{data} -\log( \mathrm{pdf}(x_\mathrm{data})) +\f] +In extended mode, a +\f$ N_\mathrm{expect} - N_\mathrm{observed}*log(N_\mathrm{expect}) \f$ term is added. +**/ + +#include + +#include "RooAbsData.h" +#include "RooAbsPdf.h" +#include "RooAbsDataStore.h" +#include "RooNLLVar.h" // RooNLLVar::ComputeScalar +#include "RunContext.h" // complete type BatchCompute::RunContext + +#include "Math/Util.h" // KahanSum + +namespace RooFit { +namespace TestStatistics { + +RooUnbinnedL::RooUnbinnedL(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended, + bool useBatchedEvaluations) + : RooAbsL(RooAbsL::ClonePdfData{pdf, data}, data->numEntries(), 1, extended), + useBatchedEvaluations_(useBatchedEvaluations) +{ +} + +RooUnbinnedL::RooUnbinnedL(const RooUnbinnedL &other) + : RooAbsL(other), apply_weight_squared(other.apply_weight_squared), _first(other._first), + useBatchedEvaluations_(other.useBatchedEvaluations_) +{ +} + +////////////////////////////////////////////////////////////////////////////////// + +/// Returns true if value was changed, false otherwise. +bool RooUnbinnedL::setApplyWeightSquared(bool flag) +{ + if (apply_weight_squared != flag) { + apply_weight_squared = flag; + return true; + } + // setValueDirty(); + return false; +} + +void RooUnbinnedL::setUseBatchedEvaluations(bool flag) { + useBatchedEvaluations_ = flag; +} + +////////////////////////////////////////////////////////////////////////////////// +/// Calculate and return likelihood on subset of data from firstEvent to lastEvent +/// processed with a step size of 'stepSize'. If this an extended likelihood and +/// and the zero event is processed the extended term is added to the return +/// likelihood. +/// +ROOT::Math::KahanSum +RooUnbinnedL::evaluatePartition(Section events, std::size_t /*components_begin*/, std::size_t /*components_end*/) +{ + // Throughout the calculation, we use Kahan's algorithm for summing to + // prevent loss of precision - this is a factor four more expensive than + // straight addition, but since evaluating the PDF is usually much more + // expensive than that, we tolerate the additional cost... + ROOT::Math::KahanSum result; + double sumWeight; + + data_->store()->recalculateCache(nullptr, events.begin(N_events_), events.end(N_events_), 1, kTRUE); + + if (useBatchedEvaluations_) { + std::tie(result, sumWeight) = RooNLLVar::computeBatchedFunc(pdf_.get(), data_.get(), evalData_, normSet_.get(), apply_weight_squared, + 1, events.begin(N_events_), events.end(N_events_)); + } else { + std::tie(result, sumWeight) = RooNLLVar::computeScalarFunc(pdf_.get(), data_.get(), normSet_.get(), apply_weight_squared, + 1, events.begin(N_events_), events.end(N_events_)); + } + + // include the extended maximum likelihood term, if requested + if (extended_) { + if (apply_weight_squared) { + + // TODO: the following should also be factored out into free/static functions like RooNLLVar::Compute* + // Calculate sum of weights-squared here for extended term + Double_t sumW2; + if (useBatchedEvaluations_) { + const RooSpan eventWeights = data_->getWeightBatch(0, N_events_); + if (eventWeights.empty()) { + sumW2 = (events.end(N_events_) - events.begin(N_events_)) * data_->weightSquared(); + } else { + ROOT::Math::KahanSum kahanWeight; + for (std::size_t i = 0; i < eventWeights.size(); ++i) { + kahanWeight.AddIndexed(eventWeights[i] * eventWeights[i], i); + } + sumW2 = kahanWeight.Sum(); + } + } else { // scalar mode + ROOT::Math::KahanSum sumW2KahanSum; + for (Int_t i = 0; i < data_->numEntries(); i++) { + data_->get(i); + sumW2KahanSum += data_->weightSquared(); + } + sumW2 = sumW2KahanSum.Sum(); + } + + Double_t expected = pdf_->expectedEvents(data_->get()); + + // Adjust calculation of extended term with W^2 weighting: adjust poisson such that + // estimate of Nexpected stays at the same value, but has a different variance, rescale + // both the observed and expected count of the Poisson with a factor sum[w] / sum[w^2] which is + // the effective weight of the Poisson term. + // i.e. change Poisson(Nobs = sum[w]| Nexp ) --> Poisson( sum[w] * sum[w] / sum[w^2] | Nexp * sum[w] / + // sum[w^2] ) weighted by the effective weight sum[w^2]/ sum[w] in the likelihood. Since here we compute + // the likelihood with the weight square we need to multiply by the square of the effective weight expectedW + // = expected * sum[w] / sum[w^2] : effective expected entries observedW = sum[w] * sum[w] / sum[w^2] : + // effective observed entries The extended term for the likelihood weighted by the square of the weight will + // be then: + // (sum[w^2]/ sum[w] )^2 * expectedW - (sum[w^2]/ sum[w] )^2 * observedW * log (expectedW) and this is + // using the previous expressions for expectedW and observedW + // sum[w^2] / sum[w] * expected - sum[w^2] * log (expectedW) + // and since the weights are constants in the likelihood we can use log(expected) instead of log(expectedW) + + Double_t expectedW2 = expected * sumW2 / data_->sumEntries(); + Double_t extra = expectedW2 - sumW2 * log(expected); + + // Double_t y = pdf->extendedTerm(sumW2, data->get()) - carry; + + result += extra; + } else { + result += pdf_->extendedTerm(data_->sumEntries(), data_->get()); + } + } + + // If part of simultaneous PDF normalize probability over + // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n) + if (sim_count_ > 1) { + result += sumWeight * log(1.0 * sim_count_); + } + + // At the end of the first full calculation, wire the caches + if (_first) { + _first = false; + pdf_->wireAllCaches(); + } + + return result; +} + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx b/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx new file mode 100644 index 0000000000000..8d61ab0697cbf --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx @@ -0,0 +1,411 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include + +#include +#include +#include +#include // necessary to complete RooUnbinnedL +#include +#include +#include +#include +#include +#include + +#include + +namespace RooFit { +/** + * \brief Namespace for new RooFit test statistic calculation. + * + * RooFit::TestStatistics contains a major refactoring of the RooAbsTestStatistic-RooAbsOptTestStatistic-RooNLLVar + * inheritance tree into: + * 1. statistics-based classes on the one hand; + * 2. calculation/evaluation/optimization based classes on the other hand. + * + * The likelihood is the central unit on the statistics side. The RooAbsL class is implemented for four kinds of + * likelihoods: binned, unbinned, "subsidiary" (an optimization for numerical stability that gathers components like + * global observables) and "sum" (over multiple components of the other types). These classes provide ways to compute + * their components in parallelizable chunks that can be used by the calculator classes as they see fit. + * + * On top of the likelihood classes, we also provide for convenience a likelihood builder buildLikelihood, as a free + * function in the namespace. This function analyzes the pdf and automatically constructs the proper likelihood, built + * up from the available RooAbsL subclasses. + * + * The calculator "Wrapper" classes are abstract interfaces. These can be implemented for different kinds of algorithms, + * or with different kinds of optimization "back-ends" in mind. In an upcoming PR, we will introduce the fork-based + * multi-processing implementation based on RooFit::MultiProcess. Other possible implementations could use the GPU or + * external tools like TensorFlow. + * + * The coupling of all these classes to RooMinimizer is made via the MinuitFcnGrad class, which owns the Wrappers that + * calculate the likelihood components. + */ +namespace TestStatistics { + +namespace { // private implementation details + +RooArgSet getConstraintsSet(RooAbsPdf *pdf, RooAbsData *data, ConstrainedParameters constrained_parameters, + ExternalConstraints external_constraints, GlobalObservables global_observables, + std::string global_observables_tag) +{ + // BEGIN CONSTRAINT COLLECTION; copied from RooAbsPdf::createNLL + + Bool_t doStripDisconnected = kFALSE; + // If no explicit list of parameters to be constrained is specified apply default algorithm + // All terms of RooProdPdfs that do not contain observables and share parameters with one or more + // terms that do contain observables are added as constraints. +#ifndef NDEBUG + bool did_default_constraint_algo = false; + std::size_t N_default_constraints = 0; +#endif + if (constrained_parameters.set.getSize() == 0) { + std::unique_ptr default_constraints{pdf->getParameters(*data, kFALSE)}; + constrained_parameters.set.add(*default_constraints); + doStripDisconnected = kTRUE; +#ifndef NDEBUG + did_default_constraint_algo = true; + N_default_constraints = default_constraints->getSize(); +#endif + } +#ifndef NDEBUG + if (did_default_constraint_algo) { + assert(N_default_constraints == static_cast(constrained_parameters.set.getSize())); + } +#endif + + // Collect internal and external constraint specifications + RooArgSet allConstraints; + + if (!global_observables_tag.empty()) { + if (global_observables.set.getSize() > 0) { + global_observables.set.removeAll(); + } + std::unique_ptr allVars{pdf->getVariables()}; + global_observables.set.add( + *dynamic_cast(allVars->selectByAttrib(global_observables_tag.c_str(), kTRUE))); + oocoutI((TObject *)nullptr, Minimization) + << "User-defined specification of global observables definition with tag named '" << global_observables_tag + << "'" << std::endl; + } else if (global_observables.set.getSize() == 0) { + // neither global_observables nor global_observables_tag was given - try if a default tag is defined in the head + // node + const char *defGlobObsTag = pdf->getStringAttribute("DefaultGlobalObservablesTag"); + if (defGlobObsTag) { + oocoutI((TObject *)nullptr, Minimization) + << "p.d.f. provides built-in specification of global observables definition with tag named '" + << defGlobObsTag << "'" << std::endl; + std::unique_ptr allVars{pdf->getVariables()}; + global_observables.set.add(*dynamic_cast(allVars->selectByAttrib(defGlobObsTag, kTRUE))); + } + } + + // EGP: removed workspace (RooAbsPdf::_myws) based stuff for now; TODO: reconnect this class to workspaces + + if (constrained_parameters.set.getSize() > 0) { + std::unique_ptr constraints{ + pdf->getAllConstraints(*data->get(), constrained_parameters.set, doStripDisconnected)}; + allConstraints.add(*constraints); + } + if (external_constraints.set.getSize() > 0) { + allConstraints.add(external_constraints.set); + } + + return allConstraints; +} + +/* + * \brief Extract a collection of subsidiary likelihoods from a pdf + * + * \param[in] pdf Raw pointer to the pdf + * \param[in] data Raw pointer to the dataset + * \param[in] constrained_parameters Set of parameters that are constrained. Pdf components dependent on these alone are + * added to the subsidiary likelihood. + * \param[in] external_constraints Set of external constraint pdfs, i.e. constraints + * not necessarily in the pdf itself. These are always added to the subsidiary likelihood. + * \param[in] global_observables + * Observables that have a constant value, independent of the dataset events. Pdf components dependent on these alone + * are added to the subsidiary likelihood. \note Overrides all other likelihood parameters (like those in \p + * constrained_parameters) if present. + * \param[in] global_observables_tag String that can be set as attribute in pdf + * components to indicate that it is a global observable. Can be used instead of or in addition to \p + * global_observables. + * \return A unique pointer to a RooSubsidiaryL that contains all terms in the pdf that can be + * calculated separately from the other components in the full likelihood. + */ +std::unique_ptr +buildSubsidiaryL(RooAbsPdf *pdf, RooAbsData *data, ConstrainedParameters constrained_parameters, + ExternalConstraints external_constraints, GlobalObservables global_observables, + std::string global_observables_tag) +{ + auto allConstraints = getConstraintsSet(pdf, data, constrained_parameters, external_constraints, global_observables, + global_observables_tag); + + std::unique_ptr subsidiary_likelihood; + // Include constraints, if any, in likelihood + if (allConstraints.getSize() > 0) { + + oocoutI((TObject *)nullptr, Minimization) + << " Including the following contraint terms in minimization: " << allConstraints << std::endl; + if (global_observables.set.getSize() > 0) { + oocoutI((TObject *)nullptr, Minimization) + << "The following global observables have been defined: " << global_observables.set << std::endl; + } + std::string name("likelihood for pdf "); + name += pdf->GetName(); + subsidiary_likelihood = std::make_unique( + name, allConstraints, + (global_observables.set.getSize() > 0) ? global_observables.set : constrained_parameters.set); + } + + return subsidiary_likelihood; +} + +bool isSimultaneous(RooAbsPdf *pdf) +{ + auto sim_pdf = dynamic_cast(pdf); + return sim_pdf != nullptr; +} + +/// Get the binned part of a pdf +/// +/// \param pdf Raw pointer to the pdf +/// \return A pdf is binned if it has attribute "BinnedLikelihood" and it is a RooRealSumPdf (or derived class). +/// If \p pdf itself is binned, it will be returned. If the pdf is a RooProdPdf (or derived), the product terms +/// will be searched for a binned component and the first such term that is found will be returned. Note that +/// the simultaneous pdf setup is such that it is assumed that only one component is binned, so this should +/// always return the correct binned component. If no binned component is found, nullptr is returned. +RooAbsPdf *getBinnedPdf(RooAbsPdf *pdf) +{ + RooAbsPdf *binnedPdf = nullptr; + if (pdf->getAttribute("BinnedLikelihood") && pdf->IsA()->InheritsFrom(RooRealSumPdf::Class())) { + // Simplest case: top-level of pdf is a RRSP + binnedPdf = pdf; + } else if (pdf->IsA()->InheritsFrom(RooProdPdf::Class())) { + // Default case: top-level pdf is a product of RRSP and other pdfs + for (const auto component : ((RooProdPdf *)pdf)->pdfList()) { + if (component->getAttribute("BinnedLikelihood") && component->IsA()->InheritsFrom(RooRealSumPdf::Class())) { + binnedPdf = (RooAbsPdf *)component; + break; + } + } + } + return binnedPdf; +} + +/* + * \brief Build a set of likelihood components to build a likelihood from a simultaneous pdf + * + * \param[in] pdf Raw pointer to the pdf + * \param[in] data Raw pointer to the dataset + * \param[in] extended Set extended term calculation on, off or use Extended::Auto to determine automatically based on + * the pdf whether to activate or not. + * \return A vector to RooAbsL unique_ptrs that contain all component binned and/or + * unbinned likelihoods. Note: subsidiary components are not included; use getConstraintsSet and/or + * buildSubsidiaryLikelihood to add those. + */ +std::vector> +getSimultaneousComponents(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended) +{ + auto sim_pdf = dynamic_cast(pdf); + + // the rest of this function is an adaptation of RooAbsTestStatistic::initSimMode: + + RooAbsCategoryLValue &simCat = (RooAbsCategoryLValue &)sim_pdf->indexCat(); + + // note: this is valid for simultaneous likelihoods, not for other test statistic types (e.g. chi2) for which this + // should return true. + bool process_empty_data_sets = RooAbsL::isExtendedHelper(pdf, extended); + + TString simCatName(simCat.GetName()); + // Note: important not to use cloned dataset here (possible when this code is run in Roo[...]L ctor), use the + // original one (which is data_ in Roo[...]L ctors, but data here) + std::unique_ptr dsetList{data->split(simCat, process_empty_data_sets)}; + if (!dsetList) { + throw std::logic_error( + "getSimultaneousComponents ERROR, index category of simultaneous pdf is missing in dataset, aborting"); + } + + // Count number of used states + std::size_t N_components = 0; + + for (const auto &catState : simCat) { + // Retrieve the PDF for this simCat state + RooAbsPdf *component_pdf = sim_pdf->getPdf(catState.first.c_str()); + auto dset = (RooAbsData *)dsetList->FindObject(catState.first.c_str()); + + if (component_pdf && dset && (0. != dset->sumEntries() || process_empty_data_sets)) { + ++N_components; + } + } + + // Allocate arrays + std::vector> components; + components.reserve(N_components); + // _gofSplitMode.resize(N_components); // not used, Hybrid mode only, see below + + // Create array of regular fit contexts, containing subset of data and single fitCat PDF + std::size_t n = 0; + for (const auto &catState : simCat) { + const std::string &catName = catState.first; + // Retrieve the PDF for this simCat state + RooAbsPdf *component_pdf = sim_pdf->getPdf(catName.c_str()); + auto dset = (RooAbsData *)dsetList->FindObject(catName.c_str()); + + if (component_pdf && dset && (0. != dset->sumEntries() || process_empty_data_sets)) { + ooccoutI((TObject *)nullptr, Fitting) + << "getSimultaneousComponents: creating slave calculator #" << n << " for state " << catName << " (" + << dset->numEntries() << " dataset entries)" << std::endl; + + RooAbsPdf *binnedPdf = getBinnedPdf(component_pdf); + Bool_t binnedL = (binnedPdf != nullptr); + if (binnedPdf == nullptr && component_pdf->IsA()->InheritsFrom(RooProdPdf::Class())) { + // Default case: top-level pdf is a product of RRSP and other pdfs + for (const auto component : ((RooProdPdf *)component_pdf)->pdfList()) { + if (component->getAttribute("MAIN_MEASUREMENT")) { + // not really a binned pdf, but this prevents a (potentially) long list of subsidiary measurements to + // be passed to the slave calculator + binnedPdf = (RooAbsPdf *)component; + break; + } + } + } + // Below here directly pass binnedPdf instead of PROD(binnedPdf,constraints) as constraints are evaluated + // elsewhere anyway and omitting them reduces model complexity and associated handling/cloning times + if (binnedL) { + components.push_back(std::make_unique((binnedPdf ? binnedPdf : component_pdf), dset)); + } else { + components.push_back(std::make_unique((binnedPdf ? binnedPdf : component_pdf), dset)); + } + // } + components.back()->setSimCount(N_components); + + // Servers may have been redirected between instantiation and (deferred) initialization + + std::unique_ptr actualParams{binnedPdf ? binnedPdf->getParameters(dset) + : component_pdf->getParameters(dset)}; + std::unique_ptr selTargetParams{ + (RooArgSet *)pdf->getParameters(*data)->selectCommon(*actualParams)}; + + assert(selTargetParams->equals(*components.back()->getParameters())); + + ++n; + } else { + if ((!dset || (0. != dset->sumEntries() && !process_empty_data_sets)) && component_pdf) { + ooccoutD((TObject *)nullptr, Fitting) << "getSimultaneousComponents: state " << catName + << " has no data entries, no slave calculator created" << std::endl; + } + } + } + oocoutI((TObject *)nullptr, Fitting) << "getSimultaneousComponents: created " << n << " slave calculators." + << std::endl; + + return components; +} + +} // anonymous namespace with private implementation details + + +/* + * \brief Build a likelihood from a pdf + dataset, optionally with a subsidiary likelihood component + * + * This function analyzes the pdf and automatically constructs the proper likelihood, built up from the available + * RooAbsL subclasses. In essence, this can give 8 conceptually different combinations, based on three questions: + * 1. Is it a simultaneous pdf? + * 2. Is the pdf binned? + * 3. Does the pdf have subsidiary terms? + * If questions 1 and 3 are answered negatively, this function will either return a RooBinnedL or RooUnbinnedL. In all + * other cases it returns a RooSumL, which will contain RooBinnedL and/or RooUnbinnedL component(s) and possibly a + * RooSubsidiaryL component with constraint terms. + * + * \param[in] pdf Raw pointer to the pdf + * \param[in] data Raw pointer to the dataset + * \param[in] extended Set extended term calculation on, off or use Extended::Auto to determine automatically based on + * the pdf whether to activate or not. + * \param[in] constrained_parameters Set of parameters that are constrained. Pdf components dependent on these alone are + * added to the subsidiary likelihood. + * \param[in] external_constraints Set of external constraint pdfs, i.e. constraints + * not necessarily in the pdf itself. These are always added to the subsidiary likelihood. + * \param[in] global_observables + * Observables that have a constant value, independent of the dataset events. Pdf components dependent on these alone + * are added to the subsidiary likelihood. \note Overrides all other likelihood parameters (like those in \p + * constrained_parameters) if present. + * \param[in] global_observables_tag String that can be set as attribute in pdf + * components to indicate that it is a global observable. Can be used instead of or in addition to \p + * global_observables. + * \return A unique pointer to a RooSubsidiaryL that contains all terms in + * the pdf that can be calculated separately from the other components in the full likelihood. + */ +std::shared_ptr buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended, + ConstrainedParameters constrained_parameters, + ExternalConstraints external_constraints, GlobalObservables global_observables, + std::string global_observables_tag) +{ + std::unique_ptr likelihood; + std::vector> components; + + if (isSimultaneous(pdf)) { + components = getSimultaneousComponents(pdf, data, extended); + } else if (auto binnedPdf = getBinnedPdf(pdf)) { + likelihood = std::make_unique(binnedPdf, data); + } else { // unbinned + likelihood = std::make_unique(pdf, data, extended); + } + + auto subsidiary = buildSubsidiaryL(pdf, data, constrained_parameters, external_constraints, global_observables, global_observables_tag); + if (subsidiary) { + if (likelihood) { + components.push_back(std::move(likelihood)); + } + components.push_back(std::move(subsidiary)); + } + if (components.size() > 0) { + likelihood = std::make_unique(pdf, data, std::move(components), extended); + } + return likelihood; +} + +// delegating convenience overloads +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, ConstrainedParameters constrained_parameters) +{ + return buildLikelihood(pdf, data, RooAbsL::Extended::Auto, constrained_parameters); +} +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, ExternalConstraints external_constraints) +{ + return buildLikelihood(pdf, data, RooAbsL::Extended::Auto, {}, external_constraints); +} +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, GlobalObservables global_observables) +{ + return buildLikelihood(pdf, data, RooAbsL::Extended::Auto, {}, {}, global_observables); +} +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, std::string global_observables_tag) +{ + return buildLikelihood(pdf, data, RooAbsL::Extended::Auto, {}, {}, {}, global_observables_tag); +} +std::shared_ptr +buildLikelihood(RooAbsPdf *pdf, RooAbsData *data, ConstrainedParameters constrained_parameters, GlobalObservables global_observables) +{ + return buildLikelihood(pdf, data, RooAbsL::Extended::Auto, constrained_parameters, {}, global_observables); +} + + +} // namespace TestStatistics +} // namespace RooFit diff --git a/roofit/roofitcore/src/TestStatistics/optional_parameter_types.cxx b/roofit/roofitcore/src/TestStatistics/optional_parameter_types.cxx new file mode 100644 index 0000000000000..599ef21fbe861 --- /dev/null +++ b/roofit/roofitcore/src/TestStatistics/optional_parameter_types.cxx @@ -0,0 +1,40 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include + +namespace RooFit { +namespace TestStatistics { + +ConstrainedParameters::ConstrainedParameters(const RooArgSet ¶meters) : set(parameters) {} +// N.B.: the default constructor must be _user-provided_ defaulted, otherwise aggregate initialization will be possible, +// which bypasses the explicit constructor and even leads to errors in some compilers; when initializing as +// ConstrainedParameters({someRooArgSet}) +// compilers can respond with +// call of overloaded ‘GlobalObservables()’ is ambiguous +// This problem is well documented in http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p1008r0.pdf. The solution +// used here is detailed in section 1.4 of that paper. Note that in C++17, the workaround is no longer necessary and +// the constructor can be _user-declared_ default (i.e. directly in the declaration above). +ConstrainedParameters::ConstrainedParameters() = default; + +ExternalConstraints::ExternalConstraints(const RooArgSet &constraints) : set(constraints) {} +ExternalConstraints::ExternalConstraints() = default; + +GlobalObservables::GlobalObservables(const RooArgSet &observables) : set(observables) {} +GlobalObservables::GlobalObservables() = default; + +} +} \ No newline at end of file diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 890723e81aefa..d64fd2252bbc8 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -44,3 +44,6 @@ ROOT_ADD_GTEST(testRooProductPdf testRooProductPdf.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore) + +ROOT_ADD_GTEST(testLikelihoodSerial TestStatistics/testLikelihoodSerial.cxx LIBRARIES RooFitCore RooFit) +ROOT_ADD_GTEST(testRooRealL TestStatistics/RooRealL.cpp LIBRARIES RooFitCore RooFit) diff --git a/roofit/roofitcore/test/TestStatistics/RooRealL.cpp b/roofit/roofitcore/test/TestStatistics/RooRealL.cpp new file mode 100644 index 0000000000000..5217d404d321d --- /dev/null +++ b/roofit/roofitcore/test/TestStatistics/RooRealL.cpp @@ -0,0 +1,336 @@ +// Author: Patrick Bos, Inti Pelupessy, Vince Croft, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include "TestStatistics/RooRealL.h" +#include "TestStatistics/RooUnbinnedL.h" +#include // necessary to complete RooUnbinnedL + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // count_if + +#include "gtest/gtest.h" + +class RooRealL + : public ::testing::TestWithParam> { +}; + +TEST_P(RooRealL, getVal) +{ + // Real-life test: calculate a NLL using event-based parallelization. This + // should replicate RooRealMPFE results. + RooRandom::randomGenerator()->SetSeed(std::get<0>(GetParam())); + RooWorkspace w; + w.factory("Gaussian::g(x[-5,5],mu[0,-3,3],sigma[1])"); + auto x = w.var("x"); + RooAbsPdf *pdf = w.pdf("g"); + std::unique_ptr data {pdf->generate(RooArgSet(*x), 10000)}; + std::unique_ptr nll {pdf->createNLL(*data)}; + + auto nominal_result = nll->getVal(); + + RooFit::TestStatistics::RooRealL nll_new("nll_new", "new style NLL", std::make_shared(pdf, data.get())); + + auto mp_result = nll_new.getVal(); + + EXPECT_DOUBLE_EQ(nominal_result, mp_result); +} + +void check_NLL_type(RooAbsReal *nll) +{ + if (dynamic_cast(nll) != nullptr) { + std::cout << "the NLL object is a RooAddition*..." << std::endl; + bool has_rooconstraintsum = false; + for (const auto nll_component : static_cast(nll)->list()) { + if (nll_component->IsA() == RooConstraintSum::Class()) { + has_rooconstraintsum = true; + std::cout << "...containing a RooConstraintSum component: " << nll_component->GetName() << std::endl; + break; + } else if (nll_component->IsA() != RooNLLVar::Class() && nll_component->IsA() != RooAddition::Class()) { + std::cerr << "... containing an unexpected component class: " << nll_component->ClassName() << std::endl; + throw std::runtime_error("RooAddition* type NLL object contains unexpected component class!"); + } + } + if (!has_rooconstraintsum) { + std::cout << "...containing only RooNLLVar components." << std::endl; + } + } else if (dynamic_cast(nll) != nullptr) { + std::cout << "the NLL object is a RooNLLVar*" << std::endl; + } +} + +void count_NLL_components(RooAbsReal *nll) +{ + if (dynamic_cast(nll) != nullptr) { + std::cout << "the NLL object is a RooAddition*..." << std::endl; + std::size_t nll_component_count = 0; + for (const auto& component : *nll->getComponents()) { + if (component->IsA() == RooNLLVar::Class()) { + ++nll_component_count; + } + } + std::cout << "...containing " << nll_component_count << " RooNLLVar components." << std::endl; + } else if (dynamic_cast(nll) != nullptr) { + std::cout << "the NLL object is a RooNLLVar*" << std::endl; + } +} + +TEST_P(RooRealL, getValRooAddition) +{ + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + RooRandom::randomGenerator()->SetSeed(std::get<0>(GetParam())); + + RooWorkspace w; + w.factory("Gaussian::g(x[-10,10],mu[0,-3,3],sigma[1])"); + + RooRealVar *x = w.var("x"); + x->setRange("x_range", -3, 0); + x->setRange("another_range", 1, 7); + + RooAbsPdf *pdf = w.pdf("g"); + std::unique_ptr data {pdf->generate(RooArgSet(*x), 10000)}; + + RooAbsReal *nll = + pdf->createNLL(*data, RooFit::Range("x_range"), RooFit::Range("another_range")); + + check_NLL_type(nll); + count_NLL_components(nll); + + delete nll; +} + +#if !defined(_MSC_VER) || defined(R__ENABLE_BROKEN_WIN_TESTS) +TEST_P(RooRealL, getValRooConstraintSumAddition) +{ + // modified from + // https://github.com/roofit-dev/rootbench/blob/43d12f33e8dac7af7d587b53a2804ddf6717e92f/root/roofit/roofit/RooFitASUM.cxx#L417 + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + int bins = 10000; + + RooRealVar x("x", "x", 0, bins); + x.setBins(bins); + // Parameters + RooRealVar a0("a0", "a0", 0); + RooRealVar a1("a1", "a1", 1, 0, 2); + RooRealVar a2("a2", "a2", 0); + + RooPolynomial p0("p0", "p0", x); + RooPolynomial p1("p1", "p1", x, RooArgList(a0, a1, a2), 0); + + std::unique_ptr dh_bkg {p0.generateBinned(x, 1000000000)}; + std::unique_ptr dh_sig {p1.generateBinned(x, 100000000)}; + dh_bkg->SetName("dh_bkg"); + dh_sig->SetName("dh_sig"); + + a1.setVal(2); + std::unique_ptr dh_sig_up {p1.generateBinned(x, 1100000000)}; + dh_sig_up->SetName("dh_sig_up"); + a1.setVal(.5); + std::unique_ptr dh_sig_down {p1.generateBinned(x, 900000000)}; + dh_sig_down->SetName("dh_sig_down"); + + RooWorkspace w = RooWorkspace("w"); + w.import(x); + w.import(*dh_sig); + w.import(*dh_bkg); + w.import(*dh_sig_up); + w.import(*dh_sig_down); + w.factory("HistFunc::hf_sig(x,dh_sig)"); + w.factory("HistFunc::hf_bkg(x,dh_bkg)"); + w.factory("HistFunc::hf_sig_up(x,dh_sig_up)"); + w.factory("HistFunc::hf_sig_down(x,dh_sig_down)"); + w.factory("PiecewiseInterpolation::pi_sig(hf_sig,hf_sig_down,hf_sig_up,alpha[-5,5])"); + + w.factory("ASUM::model(mu[1,0,5]*pi_sig,nu[1]*hf_bkg)"); + w.factory("Gaussian::constraint(alpha,0,1)"); + w.factory("PROD::model2(model,constraint)"); + + RooAbsPdf *pdf = w.pdf("model2"); + + std::unique_ptr data {pdf->generateBinned(x, 1100000)}; + RooAbsReal *nll = pdf->createNLL(*data); + + check_NLL_type(nll); + count_NLL_components(nll); + + delete nll; +} +#endif // !defined(_MSC_VER) || defined(R__ENABLE_BROKEN_WIN_TESTS) + +TEST_P(RooRealL, setVal) +{ + // calculate the NLL twice with different parameters + + RooRandom::randomGenerator()->SetSeed(std::get<0>(GetParam())); + RooWorkspace w; + w.factory("Gaussian::g(x[-5,5],mu[0,-3,3],sigma[1])"); + auto x = w.var("x"); + RooAbsPdf *pdf = w.pdf("g"); + std::unique_ptr data {pdf->generate(RooArgSet(*x), 10000)}; + std::unique_ptr nll {pdf->createNLL(*data)}; + + RooFit::TestStatistics::RooRealL nll_new("nll_new", "new style NLL", std::make_shared(pdf, data.get())); + + // calculate first results + auto nominal_result1 = nll->getVal(); + auto mp_result1 = nll_new.getVal(); + + std::cout << "nominal_result1 = " << nominal_result1 << ", mp_result1 = " << mp_result1 << std::endl; + + EXPECT_EQ(nominal_result1, mp_result1); + + w.var("mu")->setVal(2); + + // calculate second results after parameter change + auto nominal_result2 = nll->getVal(); + auto mp_result2 = nll_new.getVal(); + + std::cout << "nominal_result2 = " << nominal_result2 << ", mp_result2 = " << mp_result2 << std::endl; + + EXPECT_EQ(nominal_result2, mp_result2); + if (HasFailure()) { + std::cout << "failed test had seed = " << std::get<0>(GetParam()) << std::endl; + } +} + +INSTANTIATE_TEST_SUITE_P(NworkersModeSeed, RooRealL, + ::testing::Values(2, 3)); // random seed + +class RealLVsMPFE + : public ::testing::TestWithParam> { +}; + +TEST_P(RealLVsMPFE, getVal) +{ + // Compare our MP NLL to actual RooRealMPFE results using the same strategies. + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + // parameters + std::size_t seed = std::get<0>(GetParam()); + + RooRandom::randomGenerator()->SetSeed(seed); + + RooWorkspace w; + w.factory("Gaussian::g(x[-5,5],mu[0,-3,3],sigma[1])"); + auto x = w.var("x"); + RooAbsPdf *pdf = w.pdf("g"); + std::unique_ptr data {pdf->generate(RooArgSet(*x), 10000)}; + + std::unique_ptr nll_mpfe {pdf->createNLL(*data)}; + + auto mpfe_result = nll_mpfe->getVal(); + + RooFit::TestStatistics::RooRealL nll_new("nll_new", "new style NLL", std::make_shared(pdf, data.get())); + + auto mp_result = nll_new.getVal(); + + EXPECT_EQ(mpfe_result, mp_result); + if (HasFailure()) { + std::cout << "failed test had seed = " << std::get<0>(GetParam()) << std::endl; + } +} + +TEST_P(RealLVsMPFE, minimize) +{ + // do a minimization (e.g. like in GradMinimizer_Gaussian1D test) + + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + + // parameters + std::size_t seed = std::get<0>(GetParam()); + + RooRandom::randomGenerator()->SetSeed(seed); + + RooWorkspace w = RooWorkspace(); + + w.factory("Gaussian::g(x[-5,5],mu[0,-3,3],sigma[1])"); + auto x = w.var("x"); + RooAbsPdf *pdf = w.pdf("g"); + RooRealVar *mu = w.var("mu"); + + std::unique_ptr data {pdf->generate(RooArgSet(*x), 10000)}; + mu->setVal(-2.9); + + std::unique_ptr nll_mpfe {pdf->createNLL(*data)}; + RooFit::TestStatistics::RooRealL nll_new("nll_new", "new style NLL", std::make_shared(pdf, data.get())); + + // save initial values for the start of all minimizations + RooArgSet values = RooArgSet(*mu, *pdf); + + auto savedValues = dynamic_cast(values.snapshot()); + if (savedValues == nullptr) { + throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!"); + } + + // -------- + + RooMinimizer m0(*nll_mpfe); + m0.setMinimizerType("Minuit2"); + + m0.setStrategy(0); + m0.setPrintLevel(-1); + + m0.migrad(); + + std::unique_ptr m0result {m0.lastMinuitFit()}; + double minNll0 = m0result->minNll(); + double edm0 = m0result->edm(); + double mu0 = mu->getVal(); + double muerr0 = mu->getError(); + + values = *savedValues; + + RooMinimizer m1(nll_new); + m1.setMinimizerType("Minuit2"); + + m1.setStrategy(0); + m1.setPrintLevel(-1); + + m1.migrad(); + + std::unique_ptr m1result {m1.lastMinuitFit()}; + double minNll1 = m1result->minNll(); + double edm1 = m1result->edm(); + double mu1 = mu->getVal(); + double muerr1 = mu->getError(); + + EXPECT_EQ(minNll0, minNll1); + EXPECT_EQ(mu0, mu1); + EXPECT_EQ(muerr0, muerr1); + EXPECT_EQ(edm0, edm1); + + m1.cleanup(); // necessary in tests to clean up global _theFitter + delete savedValues; +} + +INSTANTIATE_TEST_SUITE_P(NworkersModeSeed, RealLVsMPFE, + ::testing::Values(2, 3)); // random seed diff --git a/roofit/roofitcore/test/TestStatistics/testLikelihoodSerial.cxx b/roofit/roofitcore/test/TestStatistics/testLikelihoodSerial.cxx new file mode 100644 index 0000000000000..caf5e95de1e70 --- /dev/null +++ b/roofit/roofitcore/test/TestStatistics/testLikelihoodSerial.cxx @@ -0,0 +1,443 @@ +// Author: Patrick Bos, Netherlands eScience Center / NIKHEF 2021 + +/***************************************************************************** + * RooFit + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2021, 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) * + *****************************************************************************/ + +#include + +#include +#include +#include +#include +#include +#include "RooDataHist.h" // complete type in Binned test +#include "RooCategory.h" // complete type in MultiBinnedConstraint test +#include +#include +#include +#include +#include + +#include "Math/Util.h" // KahanSum + +#include // runtime_error + +#include "gtest/gtest.h" +#include "../test_lib.h" // generate_1D_gaussian_pdf_nll + +class Environment : public testing::Environment { +public: + void SetUp() override { + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + } +}; + +// Previously, we just called AddGlobalTestEnvironment in global namespace, but this caused either a warning about an +// unused declared variable (the return value of the call) or a parsing problem that the compiler can't handle if you +// don't store the return value at all. The solution is to just define this manual main function. The default gtest +// main function does InitGoogleTest and RUN_ALL_TESTS, we add the environment call in the middle. +int main(int argc, char** argv) +{ + testing::InitGoogleTest(&argc, argv); + testing::AddGlobalTestEnvironment(new Environment); + return RUN_ALL_TESTS(); +} + +class LikelihoodSerialTest: public ::testing::Test { +protected: + void SetUp() override { + RooRandom::randomGenerator()->SetSeed(seed); + clean_flags = std::make_shared(); + } + + std::size_t seed = 23; + RooWorkspace w; + std::unique_ptr nll; + std::unique_ptr values; + RooAbsPdf *pdf; + RooAbsData *data; + std::shared_ptr likelihood; + std::shared_ptr clean_flags; +}; + +class LikelihoodSerialBinnedDatasetTest : public LikelihoodSerialTest { +protected: + void SetUp() override { + LikelihoodSerialTest::SetUp(); + + // Unbinned pdfs that define template histograms + w.factory("Gaussian::g(x[-10,10],0,2)"); + w.factory("Uniform::u(x)"); + + // Generate template histograms + RooDataHist *h_sig = w.pdf("g")->generateBinned(*w.var("x"), 1000); + RooDataHist *h_bkg = w.pdf("u")->generateBinned(*w.var("x"), 1000); + + w.import(*h_sig, RooFit::Rename("h_sig")); + w.import(*h_bkg, RooFit::Rename("h_bkg")); + + // Construct binned pdf as sum of amplitudes + w.factory("HistFunc::hf_sig(x,h_sig)"); + w.factory("HistFunc::hf_bkg(x,h_bkg)"); + w.factory("ASUM::model(mu_sig[1,-1,10]*hf_sig,mu_bkg[1,-1,10]*hf_bkg)"); + + pdf = w.pdf("model"); + } +}; + +TEST_F(LikelihoodSerialTest, UnbinnedGaussian1D) +{ + std::tie(nll, pdf, data, values) = generate_1D_gaussian_pdf_nll(w, 10000); + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialTest, UnbinnedGaussianND) +{ + unsigned int N = 4; + + std::tie(nll, pdf, data, values) = generate_ND_gaussian_pdf_nll(w, N, 1000); + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialBinnedDatasetTest, UnbinnedPdf) +{ + data = pdf->generateBinned(*w.var("x")); + + nll.reset(pdf->createNLL(*data)); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + + +TEST_F(LikelihoodSerialBinnedDatasetTest, BinnedManualNLL) +{ + pdf->setAttribute("BinnedLikelihood"); + data = pdf->generateBinned(*w.var("x")); + + // manually create NLL, ripping all relevant parts from RooAbsPdf::createNLL, except here we also set binnedL = true + RooArgSet projDeps; + RooAbsTestStatistic::Configuration nll_config; + nll_config.verbose = false; + nll_config.cloneInputData = false; + nll_config.binnedL = true; + int extended = 2; + RooNLLVar nll_manual("nlletje", "-log(likelihood)", *pdf, *data, projDeps, nll_config, extended); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll_manual.getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + + +TEST_F(LikelihoodSerialTest, SimBinned) +{ + // Unbinned pdfs that define template histograms + w.factory("Gaussian::gA(x[-10,10],-2,3)"); + w.factory("Gaussian::gB(x[-10,10],2,1)"); + w.factory("Uniform::u(x)"); + + // Generate template histograms + RooDataHist *h_sigA = w.pdf("gA")->generateBinned(*w.var("x"), 1000); + RooDataHist *h_sigB = w.pdf("gB")->generateBinned(*w.var("x"), 1000); + RooDataHist *h_bkg = w.pdf("u")->generateBinned(*w.var("x"), 1000); + + w.import(*h_sigA, RooFit::Rename("h_sigA")); + w.import(*h_sigB, RooFit::Rename("h_sigB")); + w.import(*h_bkg, RooFit::Rename("h_bkg")); + + // Construct L_phys: binned pdf as sum of amplitudes + w.factory("HistFunc::hf_sigA(x,h_sigA)"); + w.factory("HistFunc::hf_sigB(x,h_sigB)"); + w.factory("HistFunc::hf_bkg(x,h_bkg)"); + + w.factory("ASUM::model_A(mu_sig[1,-1,10]*hf_sigA,mu_bkg_A[1,-1,10]*hf_bkg)"); + w.factory("ASUM::model_B(mu_sig*hf_sigB,mu_bkg_B[1,-1,10]*hf_bkg)"); + + w.pdf("model_A")->setAttribute("BinnedLikelihood"); + w.pdf("model_B")->setAttribute("BinnedLikelihood"); + + // Construct simulatenous pdf + w.factory("SIMUL::model(index[A,B],A=model_A,B=model_B)"); + + // Construct dataset + pdf = w.pdf("model"); + data = pdf->generate(RooArgSet(*w.var("x"), *w.cat("index")), RooFit::AllBinned()); + + nll.reset(pdf->createNLL(*data)); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialTest, BinnedConstrained) +{ + // Unbinned pdfs that define template histograms + + w.factory("Gaussian::g(x[-10,10],0,2)"); + w.factory("Uniform::u(x)"); + + // Generate template histograms + + RooDataHist *h_sig = w.pdf("g")->generateBinned(*w.var("x"), 1000); + RooDataHist *h_bkg = w.pdf("u")->generateBinned(*w.var("x"), 1000); + + w.import(*h_sig, RooFit::Rename("h_sig")); + w.import(*h_bkg, RooFit::Rename("h_bkg")); + + // Construct binned pdf as sum of amplitudes + w.factory("HistFunc::hf_sig(x,h_sig)"); + w.factory("HistFunc::hf_bkg(x,h_bkg)"); + w.factory("ASUM::model_phys(mu_sig[1,-1,10]*hf_sig,expr::mu_bkg('1+0.02*alpha_bkg',alpha_bkg[-5,5])*hf_bkg)") ; + + // Construct L_subs: Gaussian subsidiary measurement that constrains alpha_bkg + w.factory("Gaussian:model_subs(alpha_bkg_obs[0],alpha_bkg,1)") ; + + // Construct full pdf + w.factory("PROD::model(model_phys,model_subs)") ; + + pdf = w.pdf("model"); + // Construct dataset from physics pdf + data = w.pdf("model_phys")->generateBinned(*w.var("x")); + + nll.reset(pdf->createNLL(*data, RooFit::GlobalObservables(*w.var("alpha_bkg_obs")))); + + // -------- + + auto nll0 = nll->getVal(); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data, RooFit::TestStatistics::GlobalObservables(RooArgSet(*w.var("alpha_bkg_obs")))); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialTest, SimUnbinned) +{ + // SIMULTANEOUS FIT OF 2 UNBINNED DATASETS + + w.factory("ExtendPdf::egA(Gaussian::gA(x[-10,10],mA[2,-10,10],s[3,0.1,10]),nA[1000])") ; + w.factory("ExtendPdf::egB(Gaussian::gB(x,mB[-2,-10,10],s),nB[100])") ; + w.factory("SIMUL::model(index[A,B],A=egA,B=egB)") ; + + pdf = w.pdf("model"); + // Construct dataset from physics pdf + data = pdf->generate(RooArgSet(*w.var("x"),*w.cat("index"))); + + nll.reset(pdf->createNLL(*data)); + + // -------- + + auto nll0 = nll->getVal(); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialTest, SimUnbinnedNonExtended) +{ + // SIMULTANEOUS FIT OF 2 UNBINNED DATASETS + + w.factory("ExtendPdf::egA(Gaussian::gA(x[-10,10],mA[2,-10,10],s[3,0.1,10]),nA[1000])") ; + w.factory("ExtendPdf::egB(Gaussian::gB(x,mB[-2,-10,10],s),nB[100])") ; + w.factory("SIMUL::model(index[A,B],A=gA,B=gB)") ; + + std::unique_ptr dA {w.pdf("gA")->generate(*w.var("x"),1)}; + std::unique_ptr dB {w.pdf("gB")->generate(*w.var("x"),1)}; + w.cat("index")->setLabel("A") ; + dA->addColumn(*w.cat("index")) ; + w.cat("index")->setLabel("B") ; + dB->addColumn(*w.cat("index")) ; + + data = (RooDataSet*) dA->Clone() ; + dynamic_cast(data)->append(*dB) ; + + pdf = w.pdf("model"); + + nll.reset(pdf->createNLL(*data)); + + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} + + +class LikelihoodSerialSimBinnedConstrainedTest : public LikelihoodSerialTest { +protected: + void SetUp() override { + LikelihoodSerialTest::SetUp(); + // Unbinned pdfs that define template histograms + + w.factory("Gaussian::gA(x[-10,10],-2,3)") ; + w.factory("Gaussian::gB(x[-10,10],2,1)") ; + w.factory("Uniform::u(x)"); + + // Generate template histograms + + std::unique_ptr h_sigA {w.pdf("gA")->generateBinned(*w.var("x"),1000)}; + std::unique_ptr h_sigB {w.pdf("gB")->generateBinned(*w.var("x"),1000)}; + std::unique_ptr h_bkg {w.pdf("u")->generateBinned(*w.var("x"), 1000)}; + + w.import(*h_sigA, RooFit::Rename("h_sigA")); + w.import(*h_sigB, RooFit::Rename("h_sigB")); + w.import(*h_bkg, RooFit::Rename("h_bkg")); + + // Construct binned pdf as sum of amplitudes + w.factory("HistFunc::hf_sigA(x,h_sigA)") ; + w.factory("HistFunc::hf_sigB(x,h_sigB)") ; + w.factory("HistFunc::hf_bkg(x,h_bkg)") ; + + w.factory("ASUM::model_phys_A(mu_sig[1,-1,10]*hf_sigA,expr::mu_bkg_A('1+0.02*alpha_bkg_A',alpha_bkg_A[-5,5])*hf_bkg)") ; + w.factory("ASUM::model_phys_B(mu_sig*hf_sigB,expr::mu_bkg_B('1+0.05*alpha_bkg_B',alpha_bkg_B[-5,5])*hf_bkg)") ; + + // Construct L_subs: Gaussian subsidiary measurement that constrains alpha_bkg + w.factory("Gaussian:model_subs_A(alpha_bkg_obs_A[0],alpha_bkg_A,1)") ; + w.factory("Gaussian:model_subs_B(alpha_bkg_obs_B[0],alpha_bkg_B,1)") ; + + // Construct full pdfs for each component (A,B) + w.factory("PROD::model_A(model_phys_A,model_subs_A)") ; + w.factory("PROD::model_B(model_phys_B,model_subs_B)") ; + + // Construct simulatenous pdf + w.factory("SIMUL::model(index[A,B],A=model_A,B=model_B)") ; + + pdf = w.pdf("model"); + // Construct dataset from physics pdf + data = pdf->generate(RooArgSet(*w.var("x"), *w.cat("index")), RooFit::AllBinned()); + } +}; + +TEST_F(LikelihoodSerialSimBinnedConstrainedTest, BasicParameters) +{ + // original test: + nll.reset(pdf->createNLL(*data, RooFit::GlobalObservables(RooArgSet(*w.var("alpha_bkg_obs_A"), *w.var("alpha_bkg_obs_B"))))); + + // -------- + + auto nll0 = nll->getVal(); + + likelihood = RooFit::TestStatistics::buildLikelihood( + pdf, data, RooFit::TestStatistics::GlobalObservables({*w.var("alpha_bkg_obs_A"), *w.var("alpha_bkg_obs_B")})); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_DOUBLE_EQ(nll0, nll1); +} + +TEST_F(LikelihoodSerialSimBinnedConstrainedTest, ConstrainedAndOffset) +{ + // a variation to test some additional parameters (ConstrainedParameters and offsetting) + nll.reset(pdf->createNLL(*data, RooFit::Constrain(RooArgSet(*w.var("alpha_bkg_obs_A"))), + RooFit::GlobalObservables(RooArgSet(*w.var("alpha_bkg_obs_B"))), RooFit::Offset(kTRUE))); + + // -------- + + auto nll0 = nll->getVal(); + + likelihood = RooFit::TestStatistics::buildLikelihood( + pdf, data, RooFit::TestStatistics::ConstrainedParameters(RooArgSet(*w.var("alpha_bkg_obs_A"))), + RooFit::TestStatistics::GlobalObservables(RooArgSet(*w.var("alpha_bkg_obs_B")))); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + nll_ts.enableOffsetting(true); + + nll_ts.evaluate(); + // The RooFit::TestStatistics classes used for minimization (RooAbsL and Wrapper derivatives) will return offset + // values, whereas RooNLLVar::getVal will always return the non-offset value, since that is the "actual" likelihood + // value. RooRealL will also give the non-offset value, so that can be directly compared to the RooNLLVar::getVal + // result (the nll0 vs nll2 comparison below). To compare to the raw RooAbsL/Wrapper value nll1, however, we need to + // manually add the offset. + ROOT::Math::KahanSum nll1 = nll_ts.getResult() + nll_ts.offset(); + + EXPECT_DOUBLE_EQ(nll0, nll1); + EXPECT_FALSE(nll_ts.offset() == 0); + + // also check against RooRealL value + RooFit::TestStatistics::RooRealL nll_real("real_nll", "RooRealL version", likelihood); + + auto nll2 = nll_real.getVal(); + + EXPECT_EQ(nll0, nll2); + EXPECT_DOUBLE_EQ(nll1, nll2); +} + +TEST_F(LikelihoodSerialTest, BatchedUnbinnedGaussianND) +{ + unsigned int N = 4; + + bool batch_mode = true; + + std::tie(nll, pdf, data, values) = generate_ND_gaussian_pdf_nll(w, N, 1000, batch_mode); + likelihood = RooFit::TestStatistics::buildLikelihood(pdf, data); + dynamic_cast(likelihood.get())->setUseBatchedEvaluations(true); + RooFit::TestStatistics::LikelihoodSerial nll_ts(likelihood, clean_flags/*, nullptr*/); + + auto nll0 = nll->getVal(); + + nll_ts.evaluate(); + auto nll1 = nll_ts.getResult(); + + EXPECT_EQ(nll0, nll1); +} diff --git a/roofit/roofitcore/test/test_lib.h b/roofit/roofitcore/test/test_lib.h index df6ead82c263d..676ade0f95642 100644 --- a/roofit/roofitcore/test/test_lib.h +++ b/roofit/roofitcore/test/test_lib.h @@ -58,7 +58,7 @@ generate_1D_gaussian_pdf_nll(RooWorkspace &w, unsigned long N_events) // return two unique_ptrs, the first because nll is a pointer, // the second because RooArgSet doesn't have a move ctor std::tuple, RooAbsPdf *, RooDataSet *, std::unique_ptr> -generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_events) { +generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_events, bool batch_mode = false) { RooArgSet obs_set; // create gaussian parameters @@ -126,7 +126,7 @@ generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_ev // --- Generate a toyMC sample from composite PDF --- RooDataSet *data = sum->generate(obs_set, N_events); - std::unique_ptr nll {sum->createNLL(*data)}; + std::unique_ptr nll {sum->createNLL(*data, RooFit::BatchMode(batch_mode))}; // set values randomly so that they actually need to do some fitting for (unsigned ix = 0; ix < n; ++ix) {