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