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