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