Logo ROOT   6.10/09
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 
35 #include "RooMsgService.h"
36 #include "TMath.h"
37 
38 #include "Riostream.h"
39 #include <assert.h>
40 
41 using 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) {
55  calculateCoefs(MaxDimension);
56  _coefsCalculated= kTRUE;
57  }
58  // allocate workspace memory
59  _nextq= new Int_t[MaxDimension];
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 {
78  _sequenceCount= 0;
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  */
101  Int_t r(0),c(_sequenceCount);
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  }
118  _sequenceCount++;
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 
205 void 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 
297 void 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 
static const Int_t _polyDegree[MaxDimension+1]
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
static const Int_t _primitivePoly[MaxDimension+1][MaxPrimitiveDegree+1]
STL namespace.
RooQuasiRandomGenerator()
Perform one-time initialization of our static coefficient array if necessary and initialize our works...
virtual ~RooQuasiRandomGenerator()
Destructor.
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Bool_t generate(UInt_t dimension, Double_t vector[])
Generate the next number in the sequence for the specified dimension.
TPaveText * pt
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
unsigned int UInt_t
Definition: RtypesCore.h:42
TMarker * m
Definition: textangle.C:8
void reset()
Reset the workspace to its initial state.
const Bool_t kFALSE
Definition: RtypesCore.h:92
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
Mother of all ROOT objects.
Definition: TObject.h:37
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.
static Int_t _cj[NBits][MaxDimension]
const Bool_t kTRUE
Definition: RtypesCore.h:91
void calculateV(const int px[], int px_degree, int pb[], int *pb_degree, int v[], int maxv)
Internal function.
void calculateCoefs(UInt_t dimension)
Calculate the coefficients for the given number of dimensions.