Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 std::endl;
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 */
95 Int_t r(0);
97 while(true) {
98 if((c % 2) == 1) {
99 ++r;
100 c /= 2;
101 }
102 else break;
103 }
104 if(r >= NBits) {
105 oocoutE(nullptr,Integration) << "RooQuasiRandomGenerator::generate: internal error!" << endl;
106 return false;
107 }
108
109 /* Calculate the next state. */
110 for(dim=0; dim<dimension; dim++) {
111 _nextq[dim] ^= _cj[r][dim];
112 }
114
115 return true;
116}
117
118
119////////////////////////////////////////////////////////////////////////////////
120/// Calculate the coefficients for the given number of dimensions
121
123{
124 int ci[NBits][NBits];
125 int v[NBits+MaxDegree+1];
126 int r;
127 unsigned int i_dim;
128
129 for(i_dim=0; i_dim<dimension; i_dim++) {
130
131 const int poly_index = i_dim + 1;
132 int j;
133 int 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; /* arbitrary 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;
216 int k;
217 int kj;
218
219 for(k=0; k<=MaxDegree; k++) {
220 ph[k] = pb[k];
221 }
222
223 /* Now multiply B by PX so B becomes PX**J.
224 * In section 2.3, the values of Bi are defined with a minus sign :
225 * don't forget this if you use them later !
226 */
227 polyMultiply(px, px_degree, pb, *pb_degree, pb, pb_degree);
228 m = *pb_degree;
229
230 /* Now choose a value of Kj as defined in section 3.3.
231 * We must have 0 <= Kj < E*J = M.
232 * The limit condition on Kj does not seem very relevant
233 * in this program.
234 */
235 /* Quoting from BFN: "Our program currently sets each K_q
236 * equal to eq. This has the effect of setting all unrestricted
237 * values of v to 1."
238 * Actually, it sets them to the arbitrary chosen value.
239 * Whatever.
240 */
241 kj = bigm;
242
243 /* Now choose values of V in accordance with
244 * the conditions in section 3.3.
245 */
246 for(r=0; r<kj; r++) {
247 v[r] = 0;
248 }
249 v[kj] = 1;
250
251
252 if(kj >= bigm) {
253 for(r=kj+1; r<m; r++) {
254 v[r] = arbitrary_element;
255 }
256 }
257 else {
258 /* This block is never reached. */
259
260 int term = sub(0, ph[kj]);
261
262 for(r=kj+1; r<bigm; r++) {
263 v[r] = arbitrary_element;
264
265 /* Check the condition of section 3.3,
266 * remembering that the H's have the opposite sign. [????????]
267 */
268 term = sub(term, mul(ph[r], v[r]));
269 }
270
271 /* Now v[bigm] != term. */
272 v[bigm] = add(nonzero_element, term);
273
274 for(r=bigm+1; r<m; r++) {
275 v[r] = arbitrary_element;
276 }
277 }
278
279 /* Calculate the remaining V's using the recursion of section 2.3,
280 * remembering that the B's have the opposite sign.
281 */
282 for(r=0; r<=maxv-m; r++) {
283 int term = 0;
284 for(k=0; k<m; k++) {
285 term = sub(term, mul(pb[k], v[r+k]));
286 }
287 v[r+m] = term;
288 }
289}
290
291
292////////////////////////////////////////////////////////////////////////////////
293/// Internal function
294
295void RooQuasiRandomGenerator::polyMultiply(const int pa[], int pa_degree, const int pb[],
296 int pb_degree, int pc[], int * pc_degree)
297{
298 int j;
299 int 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
#define c(i)
Definition RSha256.hxx:101
#define oocoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
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
TMarker m
Definition textangle.C:8