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 
22 This class generates the quasi-random (aka "low discrepancy")
23 sequence for dimensions up to 12 using the Niederreiter base 2
24 algorithm described in Bratley, Fox, Niederreiter, ACM Trans.
25 Model. Comp. Sim. 2, 195 (1992). This implementation was adapted
26 from the 0.9 beta release of the GNU scientific library.
27 Quasi-random number sequences are useful for improving the
28 convergence of a Monte Carlo integration.
29 **/
30 
31 #include "RooFit.h"
32 
34 #include "RooMsgService.h"
35 
36 #include <iostream>
37 #include <cassert>
38 
39 using 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) {
51  calculateCoefs(MaxDimension);
52  _coefsCalculated= kTRUE;
53  }
54  // allocate workspace memory
55  _nextq= new Int_t[MaxDimension];
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 {
74  _sequenceCount= 0;
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  */
97  Int_t r(0),c(_sequenceCount);
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  }
114  _sequenceCount++;
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 
201 void 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 
293 void 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 
c
#define c(i)
Definition: RSha256.hxx:119
m
auto * m
Definition: textangle.C:8
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooQuasiRandomGenerator::_polyDegree
static const Int_t _polyDegree[MaxDimension+1]
Definition: RooQuasiRandomGenerator.h:59
RooMsgService.h
RooFit.h
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
r
ROOT::R::TRInterface & r
Definition: Object.C:4
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:136
RooQuasiRandomGenerator::_primitivePoly
static const Int_t _primitivePoly[MaxDimension+1][MaxPrimitiveDegree+1]
Definition: RooQuasiRandomGenerator.h:58
oocoutE
#define oocoutE(o, a)
Definition: RooMsgService.h:48
RooQuasiRandomGenerator::MaxDimension
@ MaxDimension
Definition: RooQuasiRandomGenerator.h:41
RooQuasiRandomGenerator::RooQuasiRandomGenerator
RooQuasiRandomGenerator()
Perform one-time initialization of our static coefficient array if necessary and initialize our works...
Definition: RooQuasiRandomGenerator.cxx:48
RooQuasiRandomGenerator::~RooQuasiRandomGenerator
virtual ~RooQuasiRandomGenerator()
Destructor.
Definition: RooQuasiRandomGenerator.cxx:63
Bool_t
bool Bool_t
Definition: RtypesCore.h:63
v
@ v
Definition: rootcling_impl.cxx:3635
RooQuasiRandomGenerator::generate
Bool_t generate(UInt_t dimension, Double_t vector[])
Generate the next number in the sequence for the specified dimension.
Definition: RooQuasiRandomGenerator.cxx:83
bool
RooQuasiRandomGenerator::NBits
@ NBits
Definition: RooQuasiRandomGenerator.h:41
RooQuasiRandomGenerator
Definition: RooQuasiRandomGenerator.h:21
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
double
double
Definition: Converters.cxx:921
RooQuasiRandomGenerator::_coefsCalculated
static Bool_t _coefsCalculated
Definition: RooQuasiRandomGenerator.h:56
RooQuasiRandomGenerator::MaxPrimitiveDegree
@ MaxPrimitiveDegree
Definition: RooQuasiRandomGenerator.h:41
unsigned int
RooQuasiRandomGenerator::reset
void reset()
Reset the workspace to its initial state.
Definition: RooQuasiRandomGenerator.cxx:72
Double_t
double Double_t
Definition: RtypesCore.h:59
TObject
Definition: TObject.h:37
RooQuasiRandomGenerator::polyMultiply
void polyMultiply(const int pa[], int pa_degree, const int pb[], int pb_degree, int pc[], int *pc_degree)
Internal function.
Definition: RooQuasiRandomGenerator.cxx:293
pt
TPaveText * pt
Definition: entrylist_figure1.C:7
RooFit::Integration
@ Integration
Definition: RooGlobalFunc.h:67
RooQuasiRandomGenerator::calculateV
void calculateV(const int px[], int px_degree, int pb[], int *pb_degree, int v[], int maxv)
Internal function.
Definition: RooQuasiRandomGenerator.cxx:201
RooQuasiRandomGenerator::_cj
static Int_t _cj[NBits][MaxDimension]
Definition: RooQuasiRandomGenerator.h:57
RooQuasiRandomGenerator::calculateCoefs
void calculateCoefs(UInt_t dimension)
Calculate the coefficients for the given number of dimensions.
Definition: RooQuasiRandomGenerator.cxx:123
int
RooQuasiRandomGenerator.h