Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooAcceptReject.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooAcceptReject.cxx
19\class RooAcceptReject
20\ingroup Roofitcore
21
22Generic Monte Carlo toy generator implement
23the accept/reject sampling technique on any positively valued function.
24The RooAcceptReject generator is used by the various generator context
25classes to take care of generation of observables for which p.d.fs
26do not define internal methods
27**/
28
29#include "Riostream.h"
30
31#include "RooAcceptReject.h"
32#include "RooAbsReal.h"
33#include "RooCategory.h"
34#include "RooRealVar.h"
35#include "RooDataSet.h"
36#include "RooRandom.h"
37#include "RooErrorHandler.h"
38#include "RooPrintable.h"
39#include "RooMsgService.h"
40#include "RooRealBinding.h"
41#include "RooNumGenFactory.h"
42#include "RooNumGenConfig.h"
43
44#include "TFoam.h"
45#include "TNamed.h"
46
47#include <cassert>
48
49using std::endl, std::cerr;
50
51
52////////////////////////////////////////////////////////////////////////////////
53/// Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory
54
56{
57 RooRealVar nTrial0D("nTrial0D","Number of trial samples for cat-only generation",100,0,1e9) ;
58 RooRealVar nTrial1D("nTrial1D","Number of trial samples for 1-dim generation",1000,0,1e9) ;
59 RooRealVar nTrial2D("nTrial2D","Number of trial samples for 2-dim generation",100000,0,1e9) ;
60 RooRealVar nTrial3D("nTrial3D","Number of trial samples for N-dim generation",10000000,0,1e9) ;
61
63 fact.storeProtoSampler(proto,RooArgSet(nTrial0D,nTrial1D,nTrial2D,nTrial3D)) ;
64}
65
66
67
68////////////////////////////////////////////////////////////////////////////////
69/// Initialize an accept-reject generator for the specified distribution function,
70/// which must be non-negative but does not need to be normalized over the
71/// variables to be generated, genVars. The function and its dependents are
72/// cloned and so will not be disturbed during the generation process.
73
74RooAcceptReject::RooAcceptReject(const RooAbsReal &func, const RooArgSet &genVars, const RooNumGenConfig &config,
75 bool verbose, const RooAbsReal *maxFuncVal)
76 : RooAbsNumGenerator(func, genVars, verbose, maxFuncVal), _realSampleDim(_realVars.size()), _catSampleMult(1)
77{
78 _minTrialsArray[0] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial0D")) ;
79 _minTrialsArray[1] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial1D")) ;
80 _minTrialsArray[2] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial2D")) ;
81 _minTrialsArray[3] = static_cast<Int_t>(config.getConfigSection("RooAcceptReject").getRealValue("nTrial3D")) ;
82
83 for (auto * cat : static_range_cast<RooAbsCategory*>(_catVars)) {
84 _catSampleMult *= cat->numTypes() ;
85 }
86
87
88 // calculate the minimum number of trials needed to estimate our integral and max value
89 if (!_funcMaxVal) {
90
91 if(_realSampleDim > 3) {
93 oocoutW(nullptr, Generation) << _funcClone->GetName() << "::RooAcceptReject" << ": WARNING: generating " << _realSampleDim
94 << " variables with accept-reject may not be accurate" << endl;
95 }
96 else {
98 }
99 if (_realSampleDim > 1) {
100 oocoutW(nullptr, Generation) << "RooAcceptReject::ctor(" << _funcClone->GetName()
101 << ") WARNING: performing accept/reject sampling on a p.d.f in "
102 << _realSampleDim << " dimensions without prior knowledge on maximum value "
103 << "of p.d.f. Determining maximum value by taking " << _minTrials
104 << " trial samples. If p.d.f contains sharp peaks smaller than average "
105 << "distance between trial sampling points these may be missed and p.d.f. "
106 << "may be sampled incorrectly." << endl ;
107 }
108 } else {
109 // No trials needed if we know the maximum a priori
110 _minTrials=0 ;
111 }
112
113 // Need to fix some things here
114 if (_minTrials>10000) {
115 oocoutW(nullptr, Generation) << "RooAcceptReject::ctor(" << func.GetName() << "): WARNING: " << _minTrials << " trial samples requested by p.d.f for "
116 << _realSampleDim << "-dimensional accept/reject sampling, this may take some time" << endl ;
117 }
118
119 // print a verbose summary of our configuration, if requested
120 if(_verbose) {
121 oocoutI(nullptr, Generation) << func.GetName() << "::RooAcceptReject" << ":" << endl
122 << " Initializing accept-reject generator for" << endl << " ";
124 if (_funcMaxVal) {
125 ooccoutI(nullptr, Generation) << " Function maximum provided, no trial sampling performed" << endl ;
126 } else {
127 ooccoutI(nullptr, Generation) << " Real sampling dimension is " << _realSampleDim << endl;
128 ooccoutI(nullptr, Generation) << " Category sampling multiplier is " << _catSampleMult << endl ;
129 ooccoutI(nullptr, Generation) << " Min sampling trials is " << _minTrials << endl;
130 }
131 if (!_catVars.empty()) {
132 ooccoutI(nullptr, Generation) << " Will generate category vars "<< _catVars << endl ;
133 }
134 if (!_realVars.empty()) {
135 ooccoutI(nullptr, Generation) << " Will generate real vars " << _realVars << endl ;
136 }
137 }
138
139 // initialize our statistics
140 _maxFuncVal= 0;
141 _funcSum= 0;
142 _totalEvents= 0;
143 _eventsUsed= 0;
144}
145
146
147////////////////////////////////////////////////////////////////////////////////
148/// Return a pointer to a generated event. The caller does not own the event and it
149/// will be overwritten by a subsequent call. The input parameter 'remaining' should
150/// contain your best guess at the total number of subsequent events you will request.
151
152const RooArgSet *RooAcceptReject::generateEvent(UInt_t remaining, double& resampleRatio)
153{
154 // are we actually generating anything? (the cache always contains at least our function value)
155 const RooArgSet *event= _cache->get();
156 if(event->size() == 1) return event;
157
158 if (!_funcMaxVal) {
159 // Generation with empirical maximum determination
160
161 // first generate enough events to get reasonable estimates for the integral and
162 // maximum function value
163
164 while(_totalEvents < _minTrials) {
166
167 // Limit cache size to 1M events
168 if (_cache->numEntries()>1000000) {
169 oocoutI(nullptr, Generation) << "RooAcceptReject::generateEvent: resetting event cache" << std::endl;
170 _cache->reset() ;
171 _eventsUsed = 0 ;
172 }
173 }
174
175 event= nullptr;
176 double oldMax2(_maxFuncVal);
177 while(nullptr == event) {
178 // Use any cached events first
179 if (_maxFuncVal>oldMax2) {
180 oocxcoutD(nullptr, Generation) << "RooAcceptReject::generateEvent maxFuncVal has changed, need to resample already accepted events by factor"
181 << oldMax2 << "/" << _maxFuncVal << "=" << oldMax2/_maxFuncVal << endl ;
182 resampleRatio=oldMax2/_maxFuncVal ;
183 }
184 event= nextAcceptedEvent();
185 if(event) break;
186 // When we have used up the cache, start a new cache and add
187 // some more events to it.
188 _cache->reset();
189 _eventsUsed= 0;
190 // Calculate how many more events to generate using our best estimate of our efficiency.
191 // Always generate at least one more event so we don't get stuck.
192 if(_totalEvents*_maxFuncVal <= 0) {
193 oocoutE(nullptr, Generation) << "RooAcceptReject::generateEvent: cannot estimate efficiency...giving up" << endl;
194 return nullptr;
195 }
196
197 double eff= _funcSum/(_totalEvents*_maxFuncVal);
198 Long64_t extra= 1 + (Long64_t)(1.05*remaining/eff);
199 oocxcoutD(nullptr, Generation) << "RooAcceptReject::generateEvent: adding " << extra << " events to the cache, eff = " << eff << endl;
200 double oldMax(_maxFuncVal);
201 while(extra--) {
203 if((_maxFuncVal > oldMax)) {
204 oocxcoutD(nullptr, Generation) << "RooAcceptReject::generateEvent: estimated function maximum increased from "
205 << oldMax << " to " << _maxFuncVal << endl;
206 oldMax = _maxFuncVal ;
207 // Trim cache here
208 }
209 }
210 }
211
212 // Limit cache size to 1M events
213 if (_eventsUsed>1000000) {
214 _cache->reset() ;
215 _eventsUsed = 0 ;
216 }
217
218 } else {
219 // Generation with a priori maximum knowledge
221
222 // Generate enough trials to produce a single accepted event
223 event = nullptr ;
224 while(nullptr==event) {
226 event = nextAcceptedEvent() ;
227 }
228
229 }
230 return event;
231}
232
233
234
235////////////////////////////////////////////////////////////////////////////////
236/// Scan through events in the cache which have not been used yet,
237/// looking for the first accepted one which is added to the specified
238/// container. Return a pointer to the accepted event, or else zero
239/// if we use up the cache before we accept an event. The caller does
240/// not own the event and it will be overwritten by a subsequent call.
241
243{
244 const RooArgSet *event = nullptr;
245 while((event= _cache->get(_eventsUsed))) {
246 _eventsUsed++ ;
247 // accept this cached event?
248 double r= RooRandom::uniform();
250 //cout << " event number " << _eventsUsed << " has been rejected" << endl ;
251 continue;
252 }
253 //cout << " event number " << _eventsUsed << " has been accepted" << endl ;
254 // copy this event into the output container
255 if(_verbose && (_eventsUsed%1000==0)) {
256 cerr << "RooAcceptReject: accepted event (used " << _eventsUsed << " of "
257 << _cache->numEntries() << " so far)" << endl;
258 }
259 break;
260 }
261 //cout << "accepted event " << _eventsUsed << " of " << _cache->numEntries() << endl ;
262 return event;
263}
264
265
266
267////////////////////////////////////////////////////////////////////////////////
268/// Add a trial event to our cache and update our estimates
269/// of the function maximum value and integral.
270
272{
273 // randomize each discrete argument
274 for(auto * cat : static_range_cast<RooCategory*>(_catVars)) cat->randomize();
275
276 // randomize each real argument
277 for(auto * real : static_range_cast<RooRealVar*>(_realVars)) real->randomize();
278
279 // calculate and store our function value at this new point
280 double val= _funcClone->getVal();
281 _funcValPtr->setVal(val);
282
283 // Update the estimated integral and maximum value. Increase our
284 // maximum estimate slightly to give a safety margin with a
285 // corresponding loss of efficiency.
286 if(val > _maxFuncVal) _maxFuncVal= 1.05*val;
287 _funcSum+= val;
288
289 // fill a new entry in our cache dataset for this point
290 _cache->fill();
291 _totalEvents++;
292
293 if (_verbose &&_totalEvents%10000==0) {
294 cerr << "RooAcceptReject: generated " << _totalEvents << " events so far." << endl ;
295 }
296
297}
298
300{
301 // Empirically determine maximum value of function by taking a large number
302 // of samples. The actual number depends on the number of dimensions in which
303 // the sampling occurs
304
305 // Generate the minimum required number of samples for a reliable maximum estimate
306 while(_totalEvents < _minTrials) {
308
309 // Limit cache size to 1M events
310 if (_cache->numEntries()>1000000) {
311 oocoutI(nullptr, Generation) << "RooAcceptReject::getFuncMax: resetting event cache" << endl ;
312 _cache->reset() ;
313 _eventsUsed = 0 ;
314 }
315 }
316
317 return _maxFuncVal ;
318}
319
320std::string const& RooAcceptReject::generatorName() const {
321 static const std::string name = "RooAcceptReject";
322 return name;
323}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
#define oocxcoutD(o, a)
#define oocoutE(o, a)
#define oocoutI(o, a)
#define ooccoutI(o, a)
long long Long64_t
Definition RtypesCore.h:80
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
char name[80]
Definition TGX11.cxx:110
const char * proto
Definition civetweb.c:17536
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
Storage_t const & get() const
Const access to the underlying stl container.
Abstract base class for MC event generator implementations like RooAcceptReject and RooFoam.
RooArgSet _catVars
Set of discrete observabeles.
RooArgSet _realVars
Set of real valued observabeles.
const RooAbsReal * _funcMaxVal
Container for maximum function value.
RooAbsReal * _funcClone
Pointer to top level node of cloned function.
bool _verbose
Verbose flag.
RooRealVar * _funcValPtr
RRV storing function value in output dataset.
std::unique_ptr< RooDataSet > _cache
Dataset holding generared values of observables.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
Generic Monte Carlo toy generator implement the accept/reject sampling technique on any positively va...
UInt_t _minTrialsArray[4]
Minimum number of trials samples for 1,2,3 dimensional problems.
UInt_t _minTrials
Minimum number of max.finding trials, total number of samples.
static void registerSampler(RooNumGenFactory &fact)
Register RooIntegrator1D, is parameters and capabilities with RooNumIntFactory.
double getFuncMax() override
const RooArgSet * nextAcceptedEvent()
Scan through events in the cache which have not been used yet, looking for the first accepted one whi...
void addEventToCache()
Add a trial event to our cache and update our estimates of the function maximum value and integral.
double _funcSum
Maximum function value found, and sum of all samples made.
UInt_t _realSampleDim
Number of real dimensions to be sampled.
std::string const & generatorName() const override
Return unique name of generator implementation.
UInt_t _totalEvents
Total number of function samples.
UInt_t _catSampleMult
Number of discrete dimensions to be sampled.
const RooArgSet * generateEvent(UInt_t remaining, double &resampleRatio) override
Return a pointer to a generated event.
UInt_t _eventsUsed
Accepted number of function samples.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
Factory to instantiate numeric integrators from a given function binding and a given configuration.
bool storeProtoSampler(RooAbsNumGenerator *proto, const RooArgSet &defConfig)
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
static double uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition RooRandom.cxx:78
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47