root/RooStats/ToyMCImportanceSampler.h

0001 // @(#)root/roostats:$Id$
0002 // Author: Sven Kreiss and Kyle Cranmer    January 2012
0003 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
0004 /*************************************************************************
0005  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOSTATS_ToyMCImportanceSampler
0013 #define ROOSTATS_ToyMCImportanceSampler
0014
0015 #include "RooStats/ToyMCSampler.h"
0016 #include <vector>
0017
0018 namespace RooStats {
0019
0020 enum toysStrategies { EQUALTOYSPERDENSITY, EXPONENTIALTOYDISTRIBUTION };
0021
0022 class ToyMCImportanceSampler: public ToyMCSampler {
0023
0024    public:
0025       ToyMCImportanceSampler(TestStatistic &ts, Int_t ntoys) : ToyMCSampler(ts, ntoys) {}
0026
0027       ~ToyMCImportanceSampler() override;
0028
0029       /// overwrite GetSamplingDistributionsSingleWorker(paramPoint) with a version that loops
0030       /// over nulls and importance densities, but calls the parent
0031       /// ToyMCSampler::GetSamplingDistributionsSingleWorker(paramPoint).
0032       RooDataSet* GetSamplingDistributionsSingleWorker(RooArgSet& paramPoint) override;
0033
0034       using ToyMCSampler::GenerateToyData;
0035       RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const override;
0036       virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight, std::vector<double>& impNLLs, double& nullNLL) const;
0037       virtual RooAbsData* GenerateToyData(std::vector<double>& weights) const;
0038       virtual RooAbsData* GenerateToyData(std::vector<double>& weights, std::vector<double>& nullNLLs, std::vector<double>& impNLLs) const;
0039
0040
0041       /// specifies the pdf to sample from
0042       void SetDensityToGenerateFromByIndex(unsigned int i, bool fromNull = false) {
0043          if( (fromNull  &&  i >= fNullDensities.size())  ||
0044              (!fromNull &&  i >= fImportanceDensities.size())
0045          ) {
0046             oocoutE(nullptr,InputArguments) << "Index out of range. Requested index: "<<i<<
0047                " , but null densities: "<<fNullDensities.size()<<
0048                " and importance densities: "<<fImportanceDensities.size() << std::endl;
0049          }
0050
0051          fIndexGenDensity = i;
0052          fGenerateFromNull = fromNull;
0053
0054          ClearCache();
0055       }
0056
0057       /// For importance sampling with multiple densities/snapshots:
0058       /// This is used to check the current Likelihood against Likelihoods from
0059       /// other importance densities apart from the one given as importance snapshot.
0060       /// The pdf can be nullptr in which case the density from SetImportanceDensity()
0061       /// is used. The snapshot is also optional.
0062       void AddImportanceDensity(RooAbsPdf* p, const RooArgSet* s) {
0063          if( p == nullptr && s == nullptr ) {
0064             oocoutI(nullptr,InputArguments) << "Neither density nor snapshot given. Doing nothing." << std::endl;
0065             return;
0066          }
0067          if( p == nullptr && fPdf == nullptr ) {
0068             oocoutE(nullptr,InputArguments) << "No density given, but snapshot is there. Aborting." << std::endl;
0069             return;
0070          }
0071
0072          if( p == nullptr ) p = fPdf;
0073
0074          if( s ) s = (const RooArgSet*)s->snapshot();
0075
0076          fImportanceDensities.push_back( p );
0077          fImportanceSnapshots.push_back( s );
0078          fImpNLLs.push_back( nullptr );
0079       }
0080
0081       /// The pdf can be nullptr in which case the density from SetPdf()
0082       /// is used. The snapshot and TestStatistic is also optional.
0083       void AddNullDensity(RooAbsPdf* p, const RooArgSet* s = nullptr) {
0084          if( p == nullptr && s == nullptr ) {
0085             oocoutI(nullptr,InputArguments) << "Neither density nor snapshot nor test statistic given. Doing nothing." << std::endl;
0086             return;
0087          }
0088
0089          if( p == nullptr && !fNullDensities.empty() ) p = fNullDensities[0];
0090          if( s == nullptr ) s = fParametersForTestStat.get();
0091          if( s ) s = (const RooArgSet*)s->snapshot();
0092
0093          fNullDensities.push_back( p );
0094          fNullSnapshots.push_back( s );
0095          fNullNLLs.emplace_back( nullptr );
0096          ClearCache();
0097       }
0098       /// overwrite from ToyMCSampler
0099       void SetPdf(RooAbsPdf& pdf) override {
0100          ToyMCSampler::SetPdf(pdf);
0101
0102          if( fNullDensities.size() == 1 ) { fNullDensities[0] = &pdf; }
0103          else if( fNullDensities.empty()) AddNullDensity( &pdf );
0104          else{
0105             oocoutE(nullptr,InputArguments) << "Cannot use SetPdf() when already multiple null densities are specified. Please use AddNullDensity()." << std::endl;
0106          }
0107       }
0108       /// overwrite from ToyMCSampler
0109       void SetParametersForTestStat(const RooArgSet& nullpoi) override {
0110          ToyMCSampler::SetParametersForTestStat(nullpoi);
0111          if( fNullSnapshots.empty() ) AddNullDensity( nullptr, &nullpoi );
0112          else if( fNullSnapshots.size() == 1 ) {
0113             oocoutI(nullptr,InputArguments) << "Overwriting snapshot for the only defined null density." << std::endl;
0114             if( fNullSnapshots[0] ) delete fNullSnapshots[0];
0115             fNullSnapshots[0] = (const RooArgSet*)nullpoi.snapshot();
0116          }else{
0117             oocoutE(nullptr,InputArguments) << "Cannot use SetParametersForTestStat() when already multiple null densities are specified. Please use AddNullDensity()." << std::endl;
0118          }
0119       }
0120
0121       /// When set to true, this sets the weight of all toys to zero that
0122       /// do not have the largest likelihood under the density it was generated
0123       /// compared to the other densities.
0124       void SetApplyVeto(bool b = true) { fApplyVeto = b; }
0125
0126       void SetReuseNLL(bool r = true) { fReuseNLL = r; }
0127
0128      /// set the conditional observables which will be used when creating the NLL
0129      /// so the pdf's will not be normalized on the conditional observables when computing the NLL
0130      /// Since the class use a NLL we need to set the conditional observables if they exist in the model
0131      virtual void SetConditionalObservables(const RooArgSet& set) {fConditionalObs.removeAll(); fConditionalObs.add(set);}
0132
0133       int CreateNImpDensitiesForOnePOI(
0134          RooAbsPdf& pdf,
0135          const RooArgSet& allPOI,
0136          RooRealVar& poi,
0137          int n,
0138          double poiValueForBackground = 0.0
0139       );
0140       int CreateImpDensitiesForOnePOIAdaptively(
0141          RooAbsPdf& pdf,
0142          const RooArgSet& allPOI,
0143          RooRealVar& poi,
0144          double nStdDevOverlap = 0.5,
0145          double poiValueForBackground = 0.0
0146       );
0147
0148       void SetEqualNumToysPerDensity( void ) { fToysStrategy = EQUALTOYSPERDENSITY; }
0149       void SetExpIncreasingNumToysPerDensity( void ) { fToysStrategy = EXPONENTIALTOYDISTRIBUTION; }
0150
0151    protected:
0152
0153       /// helper method for clearing  the cache
0154       void ClearCache() override;
0155
0156       unsigned int fIndexGenDensity = 0;
0157       bool fGenerateFromNull = true;
0158       bool fApplyVeto = true;
0159
0160       RooArgSet fConditionalObs;  ///< set of conditional observables
0161
0162       /// support multiple null densities
0163       std::vector<RooAbsPdf*> fNullDensities;
0164       mutable std::vector<const RooArgSet*> fNullSnapshots;
0165
0166       // densities and snapshots to generate from
0167       std::vector<RooAbsPdf*> fImportanceDensities;
0168       std::vector<const RooArgSet*> fImportanceSnapshots;
0169
0170       bool fReuseNLL = true;
0171
0172       toysStrategies fToysStrategy = EQUALTOYSPERDENSITY;
0173
0174       mutable std::vector<std::unique_ptr<RooAbsReal>> fNullNLLs;    ///<!
0175       mutable std::vector<std::unique_ptr<RooAbsReal>> fImpNLLs;     ///<!
0176
0177    ClassDefOverride(ToyMCImportanceSampler,0) // An implementation of importance sampling
0178 };
0179 }
0180
0181 #endif