Logo ROOT  
Reference Guide
RooQuasiRandomGenerator.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 RooQuasiRandomGenerator.cxx
19\class RooQuasiRandomGenerator
20\ingroup Roofitcore
21
22This class generates the quasi-random (aka "low discrepancy")
23sequence for dimensions up to 12 using the Niederreiter base 2
24algorithm described in Bratley, Fox, Niederreiter, ACM Trans.
25Model. Comp. Sim. 2, 195 (1992). This implementation was adapted
26from the 0.9 beta release of the GNU scientific library.
27Quasi-random number sequences are useful for improving the
28convergence of a Monte Carlo integration.
29**/
30
31#include "RooFit.h"
32
34#include "RooMsgService.h"
35#include "TMath.h"
36
37#include "Riostream.h"
38#include <assert.h>
39
40using namespace std;
41
43 ;
44
45
46
47////////////////////////////////////////////////////////////////////////////////
48/// Perform one-time initialization of our static coefficient array if necessary
49/// and initialize our workspace.
50
52{
53 if(!_coefsCalculated) {
56 }
57 // allocate workspace memory
59 reset();
60}
61
62
63////////////////////////////////////////////////////////////////////////////////
64/// Destructor
65
67{
68 delete[] _nextq;
69}
70
71
72////////////////////////////////////////////////////////////////////////////////
73/// Reset the workspace to its initial state.
74
76{
78 for(Int_t dim= 0; dim < MaxDimension; dim++) _nextq[dim]= 0;
79}
80
81
82////////////////////////////////////////////////////////////////////////////////
83/// Generate the next number in the sequence for the specified dimension.
84/// The maximum dimension supported is 12.
85
87{
88 /* Load the result from the saved state. */
89 static const Double_t recip = 1.0/(double)(1U << NBits); /* 2^(-nbits) */
90
91 UInt_t dim;
92 for(dim=0; dim < dimension; dim++) {
93 vector[dim] = _nextq[dim] * recip;
94 }
95
96 /* Find the position of the least-significant zero in sequence_count.
97 * This is the bit that changes in the Gray-code representation as
98 * the count is advanced.
99 */
101 while(1) {
102 if((c % 2) == 1) {
103 ++r;
104 c /= 2;
105 }
106 else break;
107 }
108 if(r >= NBits) {
109 oocoutE((TObject*)0,Integration) << "RooQuasiRandomGenerator::generate: internal error!" << endl;
110 return kFALSE;
111 }
112
113 /* Calculate the next state. */
114 for(dim=0; dim<dimension; dim++) {
115 _nextq[dim] ^= _cj[r][dim];
116 }
118
119 return kTRUE;
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Calculate the coefficients for the given number of dimensions
125
127{
128 int ci[NBits][NBits];
129 int v[NBits+MaxDegree+1];
130 int r;
131 unsigned int i_dim;
132
133 for(i_dim=0; i_dim<dimension; i_dim++) {
134
135 const int poly_index = i_dim + 1;
136 int j, k;
137
138 /* Niederreiter (page 56, after equation (7), defines two
139 * variables Q and U. We do not need Q explicitly, but we
140 * do need U.
141 */
142 int u = 0;
143
144 /* For each dimension, we need to calculate powers of an
145 * appropriate irreducible polynomial, see Niederreiter
146 * page 65, just below equation (19).
147 * Copy the appropriate irreducible polynomial into PX,
148 * and its degree into E. Set polynomial B = PX ** 0 = 1.
149 * M is the degree of B. Subsequently B will hold higher
150 * powers of PX.
151 */
152 int pb[MaxDegree+1];
153 int px[MaxDegree+1];
154 int px_degree = _polyDegree[poly_index];
155 int pb_degree = 0;
156
157 for(k=0; k<=px_degree; k++) {
158 px[k] = _primitivePoly[poly_index][k];
159 pb[k] = 0;
160 }
161 pb[0] = 1;
162
163 for(j=0; j<NBits; j++) {
164
165 /* If U = 0, we need to set B to the next power of PX
166 * and recalculate V.
167 */
168 if(u == 0) calculateV(px, px_degree, pb, &pb_degree, v, NBits+MaxDegree);
169
170 /* Now C is obtained from V. Niederreiter
171 * obtains A from V (page 65, near the bottom), and then gets
172 * C from A (page 56, equation (7)). However this can be done
173 * in one step. Here CI(J,R) corresponds to
174 * Niederreiter's C(I,J,R).
175 */
176 for(r=0; r<NBits; r++) {
177 ci[r][j] = v[r+u];
178 }
179
180 /* Advance Niederreiter's state variables. */
181 ++u;
182 if(u == px_degree) u = 0;
183 }
184
185 /* The array CI now holds the values of C(I,J,R) for this value
186 * of I. We pack them into array CJ so that CJ(I,R) holds all
187 * the values of C(I,J,R) for J from 1 to NBITS.
188 */
189 for(r=0; r<NBits; r++) {
190 int term = 0;
191 for(j=0; j<NBits; j++) {
192 term = 2*term + ci[r][j];
193 }
194 _cj[r][i_dim] = term;
195 }
196
197 }
198}
199
200
201////////////////////////////////////////////////////////////////////////////////
202/// Internal function
203
204void RooQuasiRandomGenerator::calculateV(const int px[], int px_degree,
205 int pb[], int * pb_degree, int v[], int maxv)
206{
207 const int nonzero_element = 1; /* nonzero element of Z_2 */
208 const int arbitrary_element = 1; /* arbitray element of Z_2 */
209
210 /* The polynomial ph is px**(J-1), which is the value of B on arrival.
211 * In section 3.3, the values of Hi are defined with a minus sign:
212 * don't forget this if you use them later !
213 */
214 int ph[MaxDegree+1];
215 /* int ph_degree = *pb_degree; */
216 int bigm = *pb_degree; /* m from section 3.3 */
217 int m; /* m from section 2.3 */
218 int r, k, kj;
219
220 for(k=0; k<=MaxDegree; k++) {
221 ph[k] = pb[k];
222 }
223
224 /* Now multiply B by PX so B becomes PX**J.
225 * In section 2.3, the values of Bi are defined with a minus sign :
226 * don't forget this if you use them later !
227 */
228 polyMultiply(px, px_degree, pb, *pb_degree, pb, pb_degree);
229 m = *pb_degree;
230
231 /* Now choose a value of Kj as defined in section 3.3.
232 * We must have 0 <= Kj < E*J = M.
233 * The limit condition on Kj does not seem very relevant
234 * in this program.
235 */
236 /* Quoting from BFN: "Our program currently sets each K_q
237 * equal to eq. This has the effect of setting all unrestricted
238 * values of v to 1."
239 * Actually, it sets them to the arbitrary chosen value.
240 * Whatever.
241 */
242 kj = bigm;
243
244 /* Now choose values of V in accordance with
245 * the conditions in section 3.3.
246 */
247 for(r=0; r<kj; r++) {
248 v[r] = 0;
249 }
250 v[kj] = 1;
251
252
253 if(kj >= bigm) {
254 for(r=kj+1; r<m; r++) {
255 v[r] = arbitrary_element;
256 }
257 }
258 else {
259 /* This block is never reached. */
260
261 int term = sub(0, ph[kj]);
262
263 for(r=kj+1; r<bigm; r++) {
264 v[r] = arbitrary_element;
265
266 /* Check the condition of section 3.3,
267 * remembering that the H's have the opposite sign. [????????]
268 */
269 term = sub(term, mul(ph[r], v[r]));
270 }
271
272 /* Now v[bigm] != term. */
273 v[bigm] = add(nonzero_element, term);
274
275 for(r=bigm+1; r<m; r++) {
276 v[r] = arbitrary_element;
277 }
278 }
279
280 /* Calculate the remaining V's using the recursion of section 2.3,
281 * remembering that the B's have the opposite sign.
282 */
283 for(r=0; r<=maxv-m; r++) {
284 int term = 0;
285 for(k=0; k<m; k++) {
286 term = sub(term, mul(pb[k], v[r+k]));
287 }
288 v[r+m] = term;
289 }
290}
291
292
293////////////////////////////////////////////////////////////////////////////////
294/// Internal function
295
296void RooQuasiRandomGenerator::polyMultiply(const int pa[], int pa_degree, const int pb[],
297 int pb_degree, int pc[], int * pc_degree)
298{
299 int j, k;
300 int pt[MaxDegree+1];
301 int pt_degree = pa_degree + pb_degree;
302
303 for(k=0; k<=pt_degree; k++) {
304 int term = 0;
305 for(j=0; j<=k; j++) {
306 const int conv_term = mul(pa[k-j], pb[j]);
307 term = add(term, conv_term);
308 }
309 pt[k] = term;
310 }
311
312 for(k=0; k<=pt_degree; k++) {
313 pc[k] = pt[k];
314 }
315 for(k=pt_degree+1; k<=MaxDegree; k++) {
316 pc[k] = 0;
317 }
318
319 *pc_degree = pt_degree;
320}
321
322
323////////////////////////////////////////////////////////////////////////////////
324
327
328/* primitive polynomials in binary encoding */
329
330////////////////////////////////////////////////////////////////////////////////
331
333
334////////////////////////////////////////////////////////////////////////////////
335
337{
338 { 1, 0, 0, 0, 0, 0 }, /* 1 */
339 { 0, 1, 0, 0, 0, 0 }, /* x */
340 { 1, 1, 0, 0, 0, 0 }, /* 1 + x */
341 { 1, 1, 1, 0, 0, 0 }, /* 1 + x + x^2 */
342 { 1, 1, 0, 1, 0, 0 }, /* 1 + x + x^3 */
343 { 1, 0, 1, 1, 0, 0 }, /* 1 + x^2 + x^3 */
344 { 1, 1, 0, 0, 1, 0 }, /* 1 + x + x^4 */
345 { 1, 0, 0, 1, 1, 0 }, /* 1 + x^3 + x^4 */
346 { 1, 1, 1, 1, 1, 0 }, /* 1 + x + x^2 + x^3 + x^4 */
347 { 1, 0, 1, 0, 0, 1 }, /* 1 + x^2 + x^5 */
348 { 1, 0, 0, 1, 0, 1 }, /* 1 + x^3 + x^5 */
349 { 1, 1, 1, 1, 0, 1 }, /* 1 + x + x^2 + x^3 + x^5 */
350 { 1, 1, 1, 0, 1, 1 } /* 1 + x + x^2 + x^4 + x^5 */
351};
352
353/* degrees of primitive polynomials */
354
355////////////////////////////////////////////////////////////////////////////////
356
358{
359 0, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5
360};
361
double
Definition: Converters.cxx:921
ROOT::R::TRInterface & r
Definition: Object.C:4
#define c(i)
Definition: RSha256.hxx:101
#define oocoutE(o, a)
Definition: RooMsgService.h:48
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
This class generates the quasi-random (aka "low discrepancy") sequence for dimensions up to 12 using ...
void polyMultiply(const int pa[], int pa_degree, const int pb[], int pb_degree, int pc[], int *pc_degree)
Internal function.
RooQuasiRandomGenerator()
Perform one-time initialization of our static coefficient array if necessary and initialize our works...
void calculateV(const int px[], int px_degree, int pb[], int *pb_degree, int v[], int maxv)
Internal function.
Bool_t generate(UInt_t dimension, Double_t vector[])
Generate the next number in the sequence for the specified dimension.
Int_t mul(Int_t x, Int_t y) const
void reset()
Reset the workspace to its initial state.
void calculateCoefs(UInt_t dimension)
Calculate the coefficients for the given number of dimensions.
static Int_t _cj[NBits][MaxDimension]
virtual ~RooQuasiRandomGenerator()
Destructor.
static const Int_t _polyDegree[MaxDimension+1]
static const Int_t _primitivePoly[MaxDimension+1][MaxPrimitiveDegree+1]
Int_t add(Int_t x, Int_t y) const
Int_t sub(Int_t x, Int_t y) const
Mother of all ROOT objects.
Definition: TObject.h:37
TPaveText * pt
@ Integration
Definition: RooGlobalFunc.h:67
static constexpr double pc
auto * m
Definition: textangle.C:8