ROOT  6.06/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 //
19 // BEGIN_HTML
20 // This class generates the quasi-random (aka "low discrepancy")
21 // sequence for dimensions up to 12 using the Niederreiter base 2
22 // algorithm described in Bratley, Fox, Niederreiter, ACM Trans.
23 // Model. Comp. Sim. 2, 195 (1992). This implementation was adapted
24 // from the 0.9 beta release of the GNU scientific library.
25 // Quasi-random number sequences are useful for improving the
26 // convergence of a Monte Carlo integration.
27 // END_HTML
28 //
29 
30 #include "RooFit.h"
31 
34 #include "RooMsgService.h"
35 #include "TMath.h"
36 
37 #include "Riostream.h"
38 #include <assert.h>
39 
40 using namespace std;
41 
43  ;
44 
45 
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 /// Perform one-time initialization of our static coefficient array if necessary
49 /// and initialize our workspace.
50 
52 {
53  if(!_coefsCalculated) {
54  calculateCoefs(MaxDimension);
55  _coefsCalculated= kTRUE;
56  }
57  // allocate workspace memory
58  _nextq= new Int_t[MaxDimension];
59  reset();
60 }
61 
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 /// Destructor
65 
67 {
68  delete[] _nextq;
69 }
70 
71 
72 ////////////////////////////////////////////////////////////////////////////////
73 /// Reset the workspace to its initial state.
74 
76 {
77  _sequenceCount= 0;
78  for(Int_t dim= 0; dim < MaxDimension; dim++) _nextq[dim]= 0;
79 }
80 
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 /// Generate the next number in the sequence for the specified dimension.
84 /// The maximum dimension supported is 12.
85 
87 {
88  /* Load the result from the saved state. */
89  static const Double_t recip = 1.0/(double)(1U << NBits); /* 2^(-nbits) */
90 
91  UInt_t dim;
92  for(dim=0; dim < dimension; dim++) {
93  vector[dim] = _nextq[dim] * recip;
94  }
95 
96  /* Find the position of the least-significant zero in sequence_count.
97  * This is the bit that changes in the Gray-code representation as
98  * the count is advanced.
99  */
100  Int_t r(0),c(_sequenceCount);
101  while(1) {
102  if((c % 2) == 1) {
103  ++r;
104  c /= 2;
105  }
106  else break;
107  }
108  if(r >= NBits) {
109  oocoutE((TObject*)0,Integration) << "RooQuasiRandomGenerator::generate: internal error!" << endl;
110  return kFALSE;
111  }
112 
113  /* Calculate the next state. */
114  for(dim=0; dim<dimension; dim++) {
115  _nextq[dim] ^= _cj[r][dim];
116  }
117  _sequenceCount++;
118 
119  return kTRUE;
120 }
121 
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// Calculate the coefficients for the given number of dimensions
125 
127 {
128  int ci[NBits][NBits];
129  int v[NBits+MaxDegree+1];
130  int r;
131  unsigned int i_dim;
132 
133  for(i_dim=0; i_dim<dimension; i_dim++) {
134 
135  const int poly_index = i_dim + 1;
136  int j, k;
137 
138  /* Niederreiter (page 56, after equation (7), defines two
139  * variables Q and U. We do not need Q explicitly, but we
140  * do need U.
141  */
142  int u = 0;
143 
144  /* For each dimension, we need to calculate powers of an
145  * appropriate irreducible polynomial, see Niederreiter
146  * page 65, just below equation (19).
147  * Copy the appropriate irreducible polynomial into PX,
148  * and its degree into E. Set polynomial B = PX ** 0 = 1.
149  * M is the degree of B. Subsequently B will hold higher
150  * powers of PX.
151  */
152  int pb[MaxDegree+1];
153  int px[MaxDegree+1];
154  int px_degree = _polyDegree[poly_index];
155  int pb_degree = 0;
156 
157  for(k=0; k<=px_degree; k++) {
158  px[k] = _primitivePoly[poly_index][k];
159  pb[k] = 0;
160  }
161  pb[0] = 1;
162 
163  for(j=0; j<NBits; j++) {
164 
165  /* If U = 0, we need to set B to the next power of PX
166  * and recalculate V.
167  */
168  if(u == 0) calculateV(px, px_degree, pb, &pb_degree, v, NBits+MaxDegree);
169 
170  /* Now C is obtained from V. Niederreiter
171  * obtains A from V (page 65, near the bottom), and then gets
172  * C from A (page 56, equation (7)). However this can be done
173  * in one step. Here CI(J,R) corresponds to
174  * Niederreiter's C(I,J,R).
175  */
176  for(r=0; r<NBits; r++) {
177  ci[r][j] = v[r+u];
178  }
179 
180  /* Advance Niederreiter's state variables. */
181  ++u;
182  if(u == px_degree) u = 0;
183  }
184 
185  /* The array CI now holds the values of C(I,J,R) for this value
186  * of I. We pack them into array CJ so that CJ(I,R) holds all
187  * the values of C(I,J,R) for J from 1 to NBITS.
188  */
189  for(r=0; r<NBits; r++) {
190  int term = 0;
191  for(j=0; j<NBits; j++) {
192  term = 2*term + ci[r][j];
193  }
194  _cj[r][i_dim] = term;
195  }
196 
197  }
198 }
199 
200 
201 ////////////////////////////////////////////////////////////////////////////////
202 /// Internal function
203 
204 void RooQuasiRandomGenerator::calculateV(const int px[], int px_degree,
205  int pb[], int * pb_degree, int v[], int maxv)
206 {
207  const int nonzero_element = 1; /* nonzero element of Z_2 */
208  const int arbitrary_element = 1; /* arbitray element of Z_2 */
209 
210  /* The polynomial ph is px**(J-1), which is the value of B on arrival.
211  * In section 3.3, the values of Hi are defined with a minus sign:
212  * don't forget this if you use them later !
213  */
214  int ph[MaxDegree+1];
215  /* int ph_degree = *pb_degree; */
216  int bigm = *pb_degree; /* m from section 3.3 */
217  int m; /* m from section 2.3 */
218  int r, k, kj;
219 
220  for(k=0; k<=MaxDegree; k++) {
221  ph[k] = pb[k];
222  }
223 
224  /* Now multiply B by PX so B becomes PX**J.
225  * In section 2.3, the values of Bi are defined with a minus sign :
226  * don't forget this if you use them later !
227  */
228  polyMultiply(px, px_degree, pb, *pb_degree, pb, pb_degree);
229  m = *pb_degree;
230 
231  /* Now choose a value of Kj as defined in section 3.3.
232  * We must have 0 <= Kj < E*J = M.
233  * The limit condition on Kj does not seem very relevant
234  * in this program.
235  */
236  /* Quoting from BFN: "Our program currently sets each K_q
237  * equal to eq. This has the effect of setting all unrestricted
238  * values of v to 1."
239  * Actually, it sets them to the arbitrary chosen value.
240  * Whatever.
241  */
242  kj = bigm;
243 
244  /* Now choose values of V in accordance with
245  * the conditions in section 3.3.
246  */
247  for(r=0; r<kj; r++) {
248  v[r] = 0;
249  }
250  v[kj] = 1;
251 
252 
253  if(kj >= bigm) {
254  for(r=kj+1; r<m; r++) {
255  v[r] = arbitrary_element;
256  }
257  }
258  else {
259  /* This block is never reached. */
260 
261  int term = sub(0, ph[kj]);
262 
263  for(r=kj+1; r<bigm; r++) {
264  v[r] = arbitrary_element;
265 
266  /* Check the condition of section 3.3,
267  * remembering that the H's have the opposite sign. [????????]
268  */
269  term = sub(term, mul(ph[r], v[r]));
270  }
271 
272  /* Now v[bigm] != term. */
273  v[bigm] = add(nonzero_element, term);
274 
275  for(r=bigm+1; r<m; r++) {
276  v[r] = arbitrary_element;
277  }
278  }
279 
280  /* Calculate the remaining V's using the recursion of section 2.3,
281  * remembering that the B's have the opposite sign.
282  */
283  for(r=0; r<=maxv-m; r++) {
284  int term = 0;
285  for(k=0; k<m; k++) {
286  term = sub(term, mul(pb[k], v[r+k]));
287  }
288  v[r+m] = term;
289  }
290 }
291 
292 
293 ////////////////////////////////////////////////////////////////////////////////
294 /// Internal function
295 
296 void RooQuasiRandomGenerator::polyMultiply(const int pa[], int pa_degree, const int pb[],
297  int pb_degree, int pc[], int * pc_degree)
298 {
299  int j, 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 
static const Int_t _polyDegree[MaxDimension+1]
ClassImp(RooQuasiRandomGenerator)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
static const Int_t _primitivePoly[MaxDimension+1][MaxPrimitiveDegree+1]
const Bool_t kFALSE
Definition: Rtypes.h:92
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:48
Bool_t generate(UInt_t dimension, Double_t vector[])
Generate the next number in the sequence for the specified dimension.
TPaveText * pt
ROOT::R::TRInterface & r
Definition: Object.C:4
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.
double Double_t
Definition: RtypesCore.h:55
Mother of all ROOT objects.
Definition: TObject.h:58
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: Rtypes.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.