1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003, 2004 Ferdinando Ametrano
5 Copyright (C) 2007 Mark Joshi
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include "lowdiscrepancysequences.hpp"
22#include "utilities.hpp"
23#include <ql/math/statistics/discrepancystatistics.hpp>
24#include <ql/math/statistics/sequencestatistics.hpp>
25#include <ql/math/randomnumbers/faurersg.hpp>
26#include <ql/math/randomnumbers/haltonrsg.hpp>
27#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
28#include <ql/math/randomnumbers/seedgenerator.hpp>
29#include <ql/math/randomnumbers/primitivepolynomials.hpp>
30#include <ql/math/randomnumbers/randomizedlds.hpp>
31#include <ql/math/randomnumbers/randomsequencegenerator.hpp>
32#include <ql/math/randomnumbers/sobolrsg.hpp>
33#include <ql/utilities/dataformatters.hpp>
34#include <ql/math/randomnumbers/latticerules.hpp>
35#include <ql/math/randomnumbers/latticersg.hpp>
36
37//#define PRINT_ONLY
38#ifdef PRINT_ONLY
39#include <fstream>
40#endif
41
42using namespace QuantLib;
43using namespace boost::unit_test_framework;
44
45using std::fabs;
46
47void LowDiscrepancyTest::testSeedGenerator() {
48 BOOST_TEST_MESSAGE("Testing random-seed generator...");
49 SeedGenerator::instance().get();
50}
51
52void LowDiscrepancyTest::testPolynomialsModuloTwo() {
53
54 BOOST_TEST_MESSAGE("Testing " << PPMT_MAX_DIM <<
55 " primitive polynomials modulo two...");
56
57 const Size jj[] = {
58 1, 1, 2, 2, 6, 6, 18,
59 16, 48, 60, 176, 144, 630, 756,
60 1800, 2048, 7710, 7776, 27594, 24000, 84672,
61 120032, 356960, 276480, 1296000, 1719900, 4202496
62 };
63
64 Size i=0,j=0,n=0;
65 BigInteger polynomial=0;
66 while (n<PPMT_MAX_DIM || polynomial!=-1) {
67 if (polynomial==-1) {
68 ++i; // Increase degree index
69 j=0; // Reset index of polynomial in degree.
70 }
71 polynomial = PrimitivePolynomials[i][j];
72 if (polynomial==-1) {
73 --n;
74 if (j!=jj[i]) {
75 BOOST_ERROR("Only " << j << " polynomials in degree " << i+1
76 << " instead of " << jj[i]);
77 }
78 }
79 ++j; // Increase index of polynomial in degree i+1
80 ++n; // Increase overall polynomial counter
81 }
82
83}
84
85void LowDiscrepancyTest::testRandomizedLowDiscrepancySequence() {
86
87 BOOST_TEST_MESSAGE("Testing randomized low-discrepancy sequences up to "
88 "dimension " << PPMT_MAX_DIM << "...");
89
90 RandomizedLDS<SobolRsg, RandomSequenceGenerator<MersenneTwisterUniformRng> > rldsg(PPMT_MAX_DIM);
91 rldsg.nextSequence();
92 rldsg.lastSequence();
93 rldsg.nextRandomizer();
94
95 MersenneTwisterUniformRng t0;
96 SobolRsg t1(PPMT_MAX_DIM);
97 RandomSequenceGenerator<MersenneTwisterUniformRng> t2(PPMT_MAX_DIM);
98 RandomizedLDS<SobolRsg, RandomSequenceGenerator<MersenneTwisterUniformRng> > rldsg2(t1, t2);
99 rldsg2.nextSequence();
100 rldsg2.lastSequence();
101 rldsg2.nextRandomizer();
102
103 RandomizedLDS<SobolRsg, RandomSequenceGenerator<MersenneTwisterUniformRng> > rldsg3(t1);
104 rldsg3.nextSequence();
105 rldsg3.lastSequence();
106 rldsg3.nextRandomizer();
107
108}
109
110namespace
111{
112 void testRandomizedLatticeRule(LatticeRule::type name,
113 const std::string& nameString)
114 {
115 Size maxDim = 30;
116 Size N = 1024;
117 Size numberBatches = 32;
118
119 BOOST_TEST_MESSAGE("Testing randomized lattice sequences (" << nameString
120 << ") up to dimension " << maxDim << "...");
121
122 std::vector<Real> z;
123
124 LatticeRule::getRule(name, Z&: z, N);
125 LatticeRsg latticeGenerator(maxDim,
126 z,
127 N);
128
129 unsigned long seed = 12345678UL;
130 MersenneTwisterUniformRng rng( seed);
131
132 RandomSequenceGenerator<MersenneTwisterUniformRng> rsg(maxDim,
133 rng);
134
135 RandomizedLDS<LatticeRsg, RandomSequenceGenerator<MersenneTwisterUniformRng> > rldsg(latticeGenerator,rsg);
136
137 SequenceStatistics outerStats(maxDim);
138
139 for (Size i=0; i < numberBatches; ++i)
140 {
141 SequenceStatistics innerStats(maxDim);
142 for (Size j=0; j < N; ++j)
143 {
144 innerStats.add(sample: rldsg.nextSequence().value);
145 }
146 outerStats.add(sample: innerStats.mean());
147 rldsg.nextRandomizer();
148 }
149
150 std::vector<Real> means(outerStats.mean());
151 std::vector<Real> sds(outerStats.errorEstimate());
152
153 std::vector<Real> errorInSds(maxDim);
154
155 for (Size i=0; i < maxDim; ++i)
156 errorInSds[i] = (means[i]-0.5)/ sds[i];
157
158 Real tolerance = 4.0;
159
160 for (Size i=0; i < maxDim; ++i)
161 if (fabs(x: errorInSds[i] ) > tolerance)
162 BOOST_ERROR("Lattice generator" << nameString <<" returns a mean of " <<
163 means[i] << " with error equal to " << errorInSds[i]
164 << " standard deviations in dimension " << i);
165 }
166}
167
168
169void LowDiscrepancyTest::testRandomizedLattices()
170{
171 testRandomizedLatticeRule(name: LatticeRule::A, nameString: "A");
172 testRandomizedLatticeRule(name: LatticeRule::B, nameString: "B");
173 testRandomizedLatticeRule(name: LatticeRule::C, nameString: "C");
174 testRandomizedLatticeRule(name: LatticeRule::D, nameString: "D");
175
176}
177
178
179void LowDiscrepancyTest::testSobol() {
180
181 BOOST_TEST_MESSAGE("Testing Sobol sequences up to dimension "
182 << PPMT_MAX_DIM << "...");
183
184 std::vector<Real> point;
185 Real tolerance = 1.0e-15;
186
187 // testing max dimensionality
188 Size dimensionality = PPMT_MAX_DIM;
189 BigNatural seed = 123456;
190 SobolRsg rsg(dimensionality, seed);
191 Size points = 100, i;
192 for (i=0; i<points; i++) {
193 point = rsg.nextSequence().value;
194 if (point.size()!=dimensionality) {
195 BOOST_ERROR("Sobol sequence generator returns "
196 " a sequence of wrong dimensionality: " << point.size()
197 << " instead of " << dimensionality);
198 }
199 }
200
201 // testing homogeneity properties
202 dimensionality = 33;
203 seed = 123456;
204 rsg = SobolRsg(dimensionality, seed);
205 SequenceStatistics stat(dimensionality);
206 std::vector<Real> mean;
207 Size k = 0;
208 for (Integer j=1; j<5; j++) { // five cycle
209 points = Size(std::pow(x: 2.0, y: j))-1; // base 2
210 for (; k<points; k++) {
211 point = rsg.nextSequence().value;
212 stat.add(sample: point);
213 }
214 mean = stat.mean();
215 for (i=0; i<dimensionality; i++) {
216 Real error = std::fabs(x: mean[i]-0.5);
217 if (error > tolerance) {
218 BOOST_ERROR(io::ordinal(i+1) << " dimension: "
219 << std::fixed
220 << "mean (" << mean[i]
221 << ") at the end of the " << io::ordinal(j+1)
222 << " cycle in Sobol sequence is not " << 0.5
223 << std::scientific
224 << " (error = " << error << ")");
225 }
226 }
227 }
228
229 // testing first dimension (van der Corput sequence)
230 const Real vanderCorputSequenceModuloTwo[] = {
231 // first cycle (zero excluded)
232 0.50000,
233 // second cycle
234 0.75000, 0.25000,
235 // third cycle
236 0.37500, 0.87500, 0.62500, 0.12500,
237 // fourth cycle
238 0.18750, 0.68750, 0.93750, 0.43750, 0.31250, 0.81250, 0.56250, 0.06250,
239 // fifth cycle
240 0.09375, 0.59375, 0.84375, 0.34375, 0.46875, 0.96875, 0.71875, 0.21875,
241 0.15625, 0.65625, 0.90625, 0.40625, 0.28125, 0.78125, 0.53125, 0.03125
242 };
243
244 dimensionality = 1;
245 rsg = SobolRsg(dimensionality);
246 points = Size(std::pow(x: 2.0, y: 5))-1; // five cycles
247 for (i=0; i<points; i++) {
248 point = rsg.nextSequence().value;
249 Real error = std::fabs(x: point[0]-vanderCorputSequenceModuloTwo[i]);
250 if (error > tolerance) {
251 BOOST_ERROR(io::ordinal(i+1) << " draw ("
252 << std::fixed << point[0]
253 << ") in 1-D Sobol sequence is not in the "
254 << "van der Corput sequence modulo two: "
255 << "it should have been "
256 << vanderCorputSequenceModuloTwo[i]
257 << std::scientific
258 << " (error = " << error << ")");
259 }
260 }
261}
262
263void LowDiscrepancyTest::testFaure() {
264
265 BOOST_TEST_MESSAGE("Testing Faure sequences...");
266
267 std::vector<Real> point;
268 Real tolerance = 1.0e-15;
269
270 // testing "high" dimensionality
271 Size dimensionality = PPMT_MAX_DIM;
272 FaureRsg rsg(dimensionality);
273 Size points = 100, i;
274 for (i=0; i<points; i++) {
275 point = rsg.nextSequence().value;
276 if (point.size()!=dimensionality) {
277 BOOST_ERROR("Faure sequence generator returns "
278 " a sequence of wrong dimensionality: " << point.size()
279 << " instead of " << dimensionality);
280 }
281 }
282
283 // 1-dimension Faure (van der Corput sequence base 2)
284 const Real vanderCorputSequenceModuloTwo[] = {
285 // first cycle (zero excluded)
286 0.50000,
287 // second cycle
288 0.75000, 0.25000,
289 // third cycle
290 0.37500, 0.87500, 0.62500, 0.12500,
291 // fourth cycle
292 0.18750, 0.68750, 0.93750, 0.43750, 0.31250, 0.81250, 0.56250, 0.06250,
293 // fifth cycle
294 0.09375, 0.59375, 0.84375, 0.34375, 0.46875, 0.96875, 0.71875, 0.21875,
295 0.15625, 0.65625, 0.90625, 0.40625, 0.28125, 0.78125, 0.53125, 0.03125
296 };
297 dimensionality = 1;
298 rsg = FaureRsg(dimensionality);
299 points = Size(std::pow(x: 2.0, y: 5))-1; // five cycles
300 for (i=0; i<points; i++) {
301 point = rsg.nextSequence().value;
302 Real error = std::fabs(x: point[0]-vanderCorputSequenceModuloTwo[i]);
303 if (error > tolerance) {
304 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 1 ("
305 << std::fixed << point[0]
306 << ") in 3-D Faure sequence should have been "
307 << vanderCorputSequenceModuloTwo[i]
308 << std::scientific
309 << " (error = " << error << ")");
310 }
311 }
312
313 // 2nd dimension of the 2-dimensional Faure sequence
314 // (shuffled van der Corput sequence base 2)
315 // checked with the code provided with "Economic generation of
316 // low-discrepancy sequences with a b-ary gray code", by E. Thiemard
317 const Real FaureDimensionTwoOfTwo[] = {
318 // first cycle (zero excluded)
319 0.50000,
320 // second cycle
321 0.25000, 0.75000,
322 // third cycle
323 0.37500, 0.87500, 0.12500, 0.62500,
324 // fourth cycle
325 0.31250, 0.81250, 0.06250, 0.56250, 0.18750, 0.68750, 0.43750, 0.93750,
326 // fifth cycle
327 0.46875, 0.96875, 0.21875, 0.71875, 0.09375, 0.59375, 0.34375, 0.84375,
328 0.15625, 0.65625, 0.40625, 0.90625, 0.28125, 0.78125, 0.03125, 0.53125
329 };
330 dimensionality = 2;
331 rsg = FaureRsg(dimensionality);
332 points = Size(std::pow(x: 2.0, y: 5))-1; // five cycles
333 for (i=0; i<points; i++) {
334 point = rsg.nextSequence().value;
335 Real error = std::fabs(x: point[0]-vanderCorputSequenceModuloTwo[i]);
336 if (error > tolerance) {
337 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 1 ("
338 << std::fixed << point[0]
339 << ") in 3-D Faure sequence should have been "
340 << vanderCorputSequenceModuloTwo[i]
341 << std::scientific
342 << " (error = " << error << ")");
343 }
344 error = std::fabs(x: point[1]-FaureDimensionTwoOfTwo[i]);
345 if (error > tolerance) {
346 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 2 ("
347 << std::fixed << point[1]
348 << ") in 3-D Faure sequence should have been "
349 << FaureDimensionTwoOfTwo[i]
350 << std::scientific
351 << " (error = " << error << ")");
352 }
353 }
354
355 // 3-dimension Faure sequence (shuffled van der Corput sequence base 3)
356 // see "Monte Carlo Methods in Financial Engineering,"
357 // by Paul Glasserman, 2004 Springer Verlag, pag. 299
358 const Real FaureDimensionOneOfThree[] = {
359 // first cycle (zero excluded)
360 1.0/3, 2.0/3,
361 // second cycle
362 7.0/9, 1.0/9, 4.0/9, 5.0/9, 8.0/9, 2.0/9
363 };
364 const Real FaureDimensionTwoOfThree[] = {
365 // first cycle (zero excluded)
366 1.0/3, 2.0/3,
367 // second cycle
368 1.0/9, 4.0/9, 7.0/9, 2.0/9, 5.0/9, 8.0/9
369 };
370 const Real FaureDimensionThreeOfThree[] = {
371 // first cycle (zero excluded)
372 1.0/3, 2.0/3,
373 // second cycle
374 4.0/9, 7.0/9, 1.0/9, 8.0/9, 2.0/9, 5.0/9
375 };
376
377 dimensionality = 3;
378 rsg = FaureRsg(dimensionality);
379 points = Size(std::pow(x: 3.0, y: 2))-1; // three cycles
380 for (i=0; i<points; i++) {
381 point = rsg.nextSequence().value;
382 Real error = std::fabs(x: point[0]-FaureDimensionOneOfThree[i]);
383 if (error > tolerance) {
384 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 1 ("
385 << std::fixed << point[0]
386 << ") in 3-D Faure sequence should have been "
387 << FaureDimensionOneOfThree[i]
388 << std::scientific
389 << " (error = " << error << ")");
390 }
391 error = std::fabs(x: point[1]-FaureDimensionTwoOfThree[i]);
392 if (error > tolerance) {
393 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 2 ("
394 << std::fixed << point[1]
395 << ") in 3-D Faure sequence should have been "
396 << FaureDimensionTwoOfThree[i]
397 << std::scientific
398 << " (error = " << error << ")");
399 }
400 error = std::fabs(x: point[2]-FaureDimensionThreeOfThree[i]);
401 if (error > tolerance) {
402 BOOST_ERROR(io::ordinal(i+1) << " draw, dimension 3 ("
403 << std::fixed << point[2]
404 << ") in 3-D Faure sequence should have been "
405 << FaureDimensionThreeOfThree[i]
406 << std::scientific
407 << " (error = " << error << ")");
408 }
409 }
410}
411
412void LowDiscrepancyTest::testHalton() {
413
414 BOOST_TEST_MESSAGE("Testing Halton sequences...");
415
416 std::vector<Real> point;
417 Real tolerance = 1.0e-15;
418
419 // testing "high" dimensionality
420 Size dimensionality = PPMT_MAX_DIM;
421 HaltonRsg rsg(dimensionality, 0, false, false);
422 Size points = 100, i, k;
423 for (i=0; i<points; i++) {
424 point = rsg.nextSequence().value;
425 if (point.size()!=dimensionality) {
426 BOOST_ERROR("Halton sequence generator returns "
427 " a sequence of wrong dimensionality: " << point.size()
428 << " instead of " << dimensionality);
429 }
430 }
431
432 // testing first and second dimension (van der Corput sequence)
433 const Real vanderCorputSequenceModuloTwo[] = {
434 // first cycle (zero excluded)
435 0.50000,
436 // second cycle
437 0.25000, 0.75000,
438 // third cycle
439 0.12500, 0.62500, 0.37500, 0.87500,
440 // fourth cycle
441 0.06250, 0.56250, 0.31250, 0.81250, 0.18750, 0.68750, 0.43750, 0.93750,
442 // fifth cycle
443 0.03125, 0.53125, 0.28125, 0.78125, 0.15625, 0.65625, 0.40625, 0.90625,
444 0.09375, 0.59375, 0.34375, 0.84375, 0.21875, 0.71875, 0.46875, 0.96875,
445 };
446
447 dimensionality = 1;
448 rsg = HaltonRsg(dimensionality, 0, false, false);
449 points = Size(std::pow(x: 2.0, y: 5))-1; // five cycles
450 for (i=0; i<points; i++) {
451 point = rsg.nextSequence().value;
452 Real error = std::fabs(x: point[0]-vanderCorputSequenceModuloTwo[i]);
453 if (error > tolerance) {
454 BOOST_ERROR(io::ordinal(i+1) << " draw ("
455 << std::fixed << point[0]
456 << ") in 1-D Halton sequence is not in the "
457 << "van der Corput sequence modulo two: "
458 << "it should have been "
459 << vanderCorputSequenceModuloTwo[i]
460 << std::scientific
461 << " (error = " << error << ")");
462 }
463 }
464
465 const Real vanderCorputSequenceModuloThree[] = {
466 // first cycle (zero excluded)
467 1.0/3, 2.0/3,
468 // second cycle
469 1.0/9, 4.0/9, 7.0/9, 2.0/9, 5.0/9, 8.0/9,
470 // third cycle
471 1.0/27, 10.0/27, 19.0/27, 4.0/27, 13.0/27, 22.0/27,
472 7.0/27, 16.0/27, 25.0/27, 2.0/27, 11.0/27, 20.0/27,
473 5.0/27, 14.0/27, 23.0/27, 8.0/27, 17.0/27, 26.0/27
474 };
475
476 dimensionality = 2;
477 rsg = HaltonRsg(dimensionality, 0, false, false);
478 points = Size(std::pow(x: 3.0, y: 3))-1; // three cycles of the higher dimension
479 for (i=0; i<points; i++) {
480 point = rsg.nextSequence().value;
481 Real error = std::fabs(x: point[0]-vanderCorputSequenceModuloTwo[i]);
482 if (error > tolerance) {
483 BOOST_ERROR("First component of " << io::ordinal(i+1)
484 << " draw (" << std::fixed << point[0]
485 << ") in 2-D Halton sequence is not in the "
486 << "van der Corput sequence modulo two: "
487 << "it should have been "
488 << vanderCorputSequenceModuloTwo[i]
489 << std::scientific
490 << " (error = " << error << ")");
491 }
492 error = std::fabs(x: point[1]-vanderCorputSequenceModuloThree[i]);
493 if (error > tolerance) {
494 BOOST_ERROR("Second component of " << io::ordinal(i+1)
495 << " draw (" << std::fixed << point[1]
496 << ") in 2-D Halton sequence is not in the "
497 << "van der Corput sequence modulo three: "
498 << "it should have been "
499 << vanderCorputSequenceModuloThree[i]
500 << std::scientific
501 << " (error = " << error << ")");
502 }
503 }
504
505 // testing homogeneity properties
506 dimensionality = 33;
507 rsg = HaltonRsg(dimensionality, 0, false, false);
508 SequenceStatistics stat(dimensionality);
509 std::vector<Real> mean, stdev, variance, skewness, kurtosis;
510 k = 0;
511 Integer j;
512 for (j=1; j<5; j++) { // five cycle
513 points = Size(std::pow(x: 2.0, y: j))-1; // base 2
514 for (; k<points; k++) {
515 point = rsg.nextSequence().value;
516 stat.add(sample: point);
517 }
518 mean = stat.mean();
519 Real error = std::fabs(x: mean[0] - 0.5);
520 if (error > tolerance) {
521 BOOST_ERROR("First dimension mean (" << std::fixed << mean[0]
522 << ") at the end of the " << io::ordinal(j+1)
523 << " cycle in Halton sequence is not " << 0.5
524 << std::scientific
525 << " (error = " << error << ")");
526 }
527 }
528
529 // reset generator and gaussianstatistics
530 rsg = HaltonRsg(dimensionality, 0, false, false);
531 stat.reset(dimension: dimensionality);
532 k = 0;
533 for (j=1; j<3; j++) { // three cycle
534 points = Size(std::pow(x: 3.0, y: j))-1; // base 3
535 for (; k<points; k++) {
536 point = rsg.nextSequence().value;
537 stat.add(sample: point);
538 }
539 mean = stat.mean();
540 Real error = std::fabs(x: mean[1] - 0.5);
541 if (error > tolerance) {
542 BOOST_ERROR("Second dimension mean (" << std::fixed << mean[1]
543 << ") at the end of the " << io::ordinal(j+1)
544 << " cycle in Halton sequence is not " << 0.5
545 << std::scientific
546 << " (error = " << error << ")");
547 }
548 }
549
550}
551
552namespace {
553
554 const Real dim002Discr_Sobol[] = {
555 8.33e-004, 4.32e-004, 2.24e-004, 1.12e-004,
556 5.69e-005, 2.14e-005 // , null
557 };
558 const Real dim002DiscrMersenneTwis[] = {
559 8.84e-003, 5.42e-003, 5.23e-003, 4.47e-003,
560 4.75e-003, 3.11e-003, 2.97e-003
561 };
562 const Real dim002DiscrPlain_Halton[] = {
563 1.26e-003, 6.73e-004, 3.35e-004, 1.91e-004,
564 1.11e-004, 5.05e-005, 2.42e-005
565 };
566 const Real dim002DiscrRShiftHalton[] = {1.32e-003, 7.25e-004};
567 const Real dim002DiscrRStRShHalton[] = {1.35e-003, 9.43e-004};
568 const Real dim002DiscrRStartHalton[] = {1.08e-003, 6.40e-004};
569 const Real dim002Discr_Unit_Sobol[] = {
570 8.33e-004, 4.32e-004, 2.24e-004, 1.12e-004, 5.69e-005, 2.14e-005 // , null
571 };
572
573 const Real dim003Discr_Sobol[] = {
574 1.21e-003, 6.37e-004, 3.40e-004, 1.75e-004,
575 9.21e-005, 4.79e-005, 2.56e-005
576 };
577 const Real dim003DiscrMersenneTwis[] = {
578 7.02e-003, 4.94e-003, 4.82e-003, 4.91e-003,
579 3.33e-003, 2.80e-003, 2.62e-003
580 };
581 const Real dim003DiscrPlain_Halton[] = {
582 1.63e-003, 9.62e-004, 4.83e-004, 2.67e-004,
583 1.41e-004, 7.64e-005, 3.93e-005
584 };
585 const Real dim003DiscrRShiftHalton[] = {1.96e-003, 1.03e-003};
586 const Real dim003DiscrRStRShHalton[] = {2.17e-003, 1.54e-003};
587 const Real dim003DiscrRStartHalton[] = {1.48e-003, 7.77e-004};
588 const Real dim003Discr_Unit_Sobol[] = {1.21e-003, 6.37e-004, 3.40e-004, 1.75e-004,
589 9.21e-005, 4.79e-005, 2.56e-005};
590
591 const Real dim005Discr_Sobol[] = {
592 1.59e-003, 9.55e-004, 5.33e-004, 3.22e-004,
593 1.63e-004, 9.41e-005, 5.19e-005
594 };
595 const Real dim005DiscrMersenneTwis[] = {
596 4.28e-003, 3.48e-003, 2.48e-003, 1.98e-003,
597 1.57e-003, 1.39e-003, 6.33e-004
598 };
599 const Real dim005DiscrPlain_Halton[] = {
600 1.93e-003, 1.23e-003, 6.89e-004, 4.22e-004,
601 2.13e-004, 1.25e-004, 7.17e-005
602 };
603 const Real dim005DiscrRShiftHalton[] = {2.02e-003, 1.36e-003};
604 const Real dim005DiscrRStRShHalton[] = {2.11e-003, 1.25e-003};
605 const Real dim005DiscrRStartHalton[] = {1.74e-003, 1.08e-003};
606 const Real dim005Discr_Unit_Sobol[] = {1.85e-003, 9.39e-004, 5.19e-004, 2.99e-004,
607 1.75e-004, 9.51e-005, 5.55e-005};
608
609 const Real dim010DiscrJackel_Sobol[] = {
610 7.08e-004, 5.31e-004, 3.60e-004, 2.18e-004,
611 1.57e-004, 1.12e-004, 6.39e-005
612 };
613 const Real dim010DiscrSobLev_Sobol[] = {
614 7.01e-004, 5.10e-004, 3.28e-004, 2.21e-004,
615 1.57e-004, 1.08e-004, 6.38e-005
616 };
617 const Real dim010DiscrMersenneTwis[] = {
618 8.83e-004, 6.56e-004, 4.87e-004, 3.37e-004,
619 3.06e-004, 1.73e-004, 1.43e-004
620 };
621 const Real dim010DiscrPlain_Halton[] = {
622 1.23e-003, 6.89e-004, 4.03e-004, 2.83e-004,
623 1.61e-004, 1.08e-004, 6.69e-005
624 };
625 const Real dim010DiscrRShiftHalton[] = {9.25e-004, 6.40e-004};
626 const Real dim010DiscrRStRShHalton[] = {8.41e-004, 5.42e-004};
627 const Real dim010DiscrRStartHalton[] = {7.89e-004, 5.33e-004};
628 const Real dim010Discr_Unit_Sobol[] = {7.67e-004, 4.92e-004, 3.47e-004, 2.34e-004,
629 1.39e-004, 9.47e-005, 5.72e-005};
630
631 const Real dim015DiscrJackel_Sobol[] = {
632 1.59e-004, 1.23e-004, 7.73e-005, 5.51e-005,
633 3.91e-005, 2.73e-005, 1.96e-005
634 };
635 const Real dim015DiscrSobLev_Sobol[] = {
636 1.48e-004, 1.06e-004, 8.19e-005, 6.29e-005,
637 4.16e-005, 2.54e-005, 1.73e-005
638 };
639 const Real dim015DiscrMersenneTwis[] = {
640 1.63e-004, 1.12e-004, 8.36e-005, 6.09e-005,
641 4.34e-005, 2.95e-005, 2.10e-005
642 };
643 const Real dim015DiscrPlain_Halton[] = {
644 5.75e-004, 3.12e-004, 1.70e-004, 9.89e-005,
645 5.33e-005, 3.45e-005, 2.11e-005
646 };
647 const Real dim015DiscrRShiftHalton[] = {1.75e-004, 1.19e-004};
648 const Real dim015DiscrRStRShHalton[] = {1.66e-004, 1.34e-004};
649 const Real dim015DiscrRStartHalton[] = {2.09e-004, 1.30e-004};
650 const Real dim015Discr_Unit_Sobol[] = {2.24e-004, 1.39e-004, 9.86e-005, 6.02e-005,
651 4.39e-005, 3.06e-005, 2.32e-005};
652
653 const Real dim030DiscrJackel_Sobol[] = {
654 6.43e-007, 5.28e-007, 3.88e-007, 2.49e-007,
655 2.09e-007, 1.55e-007, 1.07e-007
656 };
657 const Real dim030DiscrSobLev_Sobol[] = {
658 1.03e-006, 6.06e-007, 3.81e-007, 2.71e-007,
659 2.68e-007, 1.73e-007, 1.21e-007
660 };
661 const Real dim030DiscrMersenneTwis[] = {
662 4.38e-007, 3.25e-007, 4.47e-007, 2.85e-007,
663 2.03e-007, 1.50e-007, 1.17e-007
664 };
665 const Real dim030DiscrPlain_Halton[] = {
666 4.45e-004, 2.23e-004, 1.11e-004, 5.56e-005,
667 2.78e-005, 1.39e-005, 6.95e-006
668 };
669 const Real dim030DiscrRShiftHalton[] = {8.11e-007, 6.05e-007};
670 const Real dim030DiscrRStRShHalton[] = {1.85e-006, 1.03e-006};
671 const Real dim030DiscrRStartHalton[] = {4.42e-007, 4.64e-007};
672 const Real dim030Discr_Unit_Sobol[] = {4.35e-005, 2.17e-005, 1.09e-005, 5.43e-006,
673 2.73e-006, 1.37e-006, 6.90e-007};
674
675 const Real dim050DiscrJackel_Sobol[] = {
676 2.98e-010, 2.91e-010, 2.62e-010, 1.53e-010,
677 1.48e-010, 1.15e-010, 8.41e-011
678 };
679 const Real dim050DiscrSobLev_Sobol[] = {
680 3.11e-010, 2.52e-010, 1.61e-010, 1.54e-010,
681 1.11e-010, 8.60e-011, 1.17e-010
682 };
683 const Real dim050DiscrSobLem_Sobol[] = {
684 4.57e-010, 6.84e-010, 3.68e-010, 2.20e-010,
685 1.81e-010, 1.14e-010, 8.31e-011
686 };
687 const Real dim050DiscrMersenneTwis[] = {
688 3.27e-010, 2.42e-010, 1.47e-010, 1.98e-010,
689 2.31e-010, 1.30e-010, 8.09e-011
690 };
691 const Real dim050DiscrPlain_Halton[] = {
692 4.04e-004, 2.02e-004, 1.01e-004, 5.05e-005,
693 2.52e-005, 1.26e-005, 6.31e-006
694 };
695 const Real dim050DiscrRShiftHalton[] = {1.14e-010, 1.25e-010};
696 const Real dim050DiscrRStRShHalton[] = {2.92e-010, 5.02e-010};
697 const Real dim050DiscrRStartHalton[] = {1.93e-010, 6.82e-010};
698 const Real dim050Discr_Unit_Sobol[] = {1.63e-005, 8.14e-006, 4.07e-006, 2.04e-006,
699 1.02e-006, 5.09e-007, 2.54e-007};
700
701 const Real dim100DiscrJackel_Sobol[] = {
702 1.26e-018, 1.55e-018, 8.46e-019, 4.43e-019,
703 4.04e-019, 2.44e-019, 4.86e-019
704 };
705 const Real dim100DiscrSobLev_Sobol[] = {
706 1.17e-018, 2.65e-018, 1.45e-018, 7.28e-019,
707 6.33e-019, 3.36e-019, 3.43e-019
708 };
709 const Real dim100DiscrSobLem_Sobol[] = {
710 8.79e-019, 4.60e-019, 6.69e-019, 7.17e-019,
711 5.81e-019, 2.97e-019, 2.64e-019
712 };
713 const Real dim100DiscrMersenneTwis[] = {
714 5.30e-019, 7.29e-019, 3.71e-019, 3.33e-019,
715 1.33e-017, 6.70e-018, 3.36e-018
716 };
717 const Real dim100DiscrPlain_Halton[] = {
718 3.63e-004, 1.81e-004, 9.07e-005, 4.53e-005,
719 2.27e-005, 1.13e-005, 5.66e-006
720 };
721 const Real dim100DiscrRShiftHalton[] = {3.36e-019, 2.19e-019};
722 const Real dim100DiscrRStRShHalton[] = {4.44e-019, 2.24e-019};
723 const Real dim100DiscrRStartHalton[] = {9.85e-020, 8.34e-019};
724 const Real dim100Discr_Unit_Sobol[] = {4.97e-006, 2.48e-006, 1.24e-006, 6.20e-007,
725 3.10e-007, 1.55e-007, 7.76e-008};
726
727 const Size dimensionality[] = {2, 3, 5, 10, 15, 30, 50, 100 };
728
729 // 7 discrepancy measures for each dimension of all sequence generators
730 // would take a few days ... too long for usual/frequent test running
731 const Size discrepancyMeasuresNumber = 1;
732
733 // let's add some generality here...
734
735 class MersenneFactory {
736 public:
737 typedef RandomSequenceGenerator<MersenneTwisterUniformRng>
738 MersenneTwisterUniformRsg;
739 typedef MersenneTwisterUniformRsg generator_type;
740 MersenneTwisterUniformRsg make(Size dim,
741 BigNatural seed) const {
742 return MersenneTwisterUniformRsg(dim,seed);
743 }
744 std::string name() const { return "Mersenne Twister"; }
745 };
746
747 class SobolFactory {
748 public:
749 typedef SobolRsg generator_type;
750 explicit SobolFactory(SobolRsg::DirectionIntegers unit) : unit_(unit) {}
751 SobolRsg make(Size dim,
752 BigNatural seed) const {
753 return SobolRsg(dim,seed,unit_);
754 }
755 std::string name() const {
756 std::string prefix;
757 switch (unit_) {
758 case SobolRsg::Unit:
759 prefix = "unit-initialized ";
760 break;
761 case SobolRsg::Jaeckel:
762 prefix = "Jaeckel-initialized ";
763 break;
764 case SobolRsg::SobolLevitan:
765 prefix = "SobolLevitan-initialized ";
766 break;
767 case SobolRsg::SobolLevitanLemieux:
768 prefix = "SobolLevitanLemieux-initialized ";
769 break;
770 case SobolRsg::Kuo:
771 prefix = "Kuo";
772 break;
773 case SobolRsg::Kuo2:
774 prefix = "Kuo2";
775 break;
776 case SobolRsg::Kuo3:
777 prefix = "Kuo3";
778 break;
779 default:
780 QL_FAIL("unknown direction integers");
781 }
782 return prefix + "Sobol sequences: ";
783 }
784 private:
785 SobolRsg::DirectionIntegers unit_;
786 };
787
788 class HaltonFactory {
789 public:
790 typedef HaltonRsg generator_type;
791 HaltonFactory(bool randomStart, bool randomShift)
792 : start_(randomStart), shift_(randomShift) {}
793 HaltonRsg make(Size dim,
794 BigNatural seed) const {
795 return HaltonRsg(dim,seed,start_,shift_);
796 }
797 std::string name() const {
798 std::string prefix = start_ ?
799 "random-start " :
800 "";
801 if (shift_)
802 prefix += "random-shift ";
803 return prefix + "Halton";
804 }
805 private:
806 bool start_, shift_;
807 };
808
809 template <class T>
810 void testGeneratorDiscrepancy(const T& generatorFactory,
811 #ifndef PRINT_ONLY
812 const Real * const discrepancy[8],
813 const std::string&,
814 const std::string&
815 #else
816 const Real * const [8],
817 const std::string& fileName,
818 const std::string& arrayName
819 #endif
820 ) {
821
822 #ifndef PRINT_ONLY
823 Real tolerance = 1.0e-2;
824 #endif
825
826 std::vector<Real> point;
827 Size dim;
828 BigNatural seed = 123456;
829 Real discr;
830 // more than 1 discrepancy measures take long time
831 Size sampleLoops = std::max<Size>(a: 1, b: discrepancyMeasuresNumber);
832
833 #ifdef PRINT_ONLY
834 std::ofstream outStream(fileName.c_str());
835 #endif
836 for (Integer i = 0; i<8; i++) {
837 #ifdef PRINT_ONLY
838 outStream << std::endl;
839 #endif
840
841 dim = dimensionality[i];
842 DiscrepancyStatistics stat(dim);
843
844 typename T::generator_type rsg = generatorFactory.make(dim, seed);
845
846 Size j, k=0, jMin=10;
847 stat.reset();
848 #ifdef PRINT_ONLY
849 outStream << "const Real dim" << dim
850 << arrayName << "[] = {" ;
851 #endif
852 for (j=jMin; j<jMin+sampleLoops; j++) {
853 Size points = Size(std::pow(x: 2.0, y: Integer(j)))-1;
854 for (; k<points; k++) {
855 point = rsg.nextSequence().value;
856 stat.add(sample: point);
857 }
858
859 discr = stat.discrepancy();
860
861 #ifdef PRINT_ONLY
862 if (j!=jMin)
863 outStream << ", ";
864 outStream << std::fixed << std::setprecision(2) << discr;
865 #else
866 if (std::fabs(x: discr-discrepancy[i][j-jMin])>tolerance*discr) {
867 BOOST_ERROR(generatorFactory.name()
868 << "discrepancy dimension " << dimensionality[i]
869 << " at " << points << " samples is "
870 << discr << " instead of "
871 << discrepancy[i][j-jMin]);
872 }
873 #endif
874 }
875 #ifdef PRINT_ONLY
876 outStream << "};" << std::endl;
877 #endif
878 }
879 #ifdef PRINT_ONLY
880 outStream.close();
881 #endif
882 }
883
884}
885
886
887void LowDiscrepancyTest::testMersenneTwisterDiscrepancy() {
888
889 BOOST_TEST_MESSAGE("Testing Mersenne-twister discrepancy...");
890
891 const Real * const discrepancy[8] = {
892 dim002DiscrMersenneTwis, dim003DiscrMersenneTwis,
893 dim005DiscrMersenneTwis, dim010DiscrMersenneTwis,
894 dim015DiscrMersenneTwis, dim030DiscrMersenneTwis,
895 dim050DiscrMersenneTwis, dim100DiscrMersenneTwis
896 };
897
898 testGeneratorDiscrepancy(generatorFactory: MersenneFactory(),
899 discrepancy,
900 "MersenneDiscrepancy.txt",
901 "DiscrMersenneTwis");
902}
903
904void LowDiscrepancyTest::testPlainHaltonDiscrepancy() {
905
906 BOOST_TEST_MESSAGE("Testing plain Halton discrepancy...");
907
908 const Real * const discrepancy[8] = {
909 dim002DiscrPlain_Halton, dim003DiscrPlain_Halton,
910 dim005DiscrPlain_Halton, dim010DiscrPlain_Halton,
911 dim015DiscrPlain_Halton, dim030DiscrPlain_Halton,
912 dim050DiscrPlain_Halton, dim100DiscrPlain_Halton};
913
914 testGeneratorDiscrepancy(generatorFactory: HaltonFactory(false,false),
915 discrepancy,
916 "PlainHaltonDiscrepancy.txt",
917 "DiscrPlain_Halton");
918}
919
920void LowDiscrepancyTest::testRandomStartHaltonDiscrepancy() {
921
922 BOOST_TEST_MESSAGE("Testing random-start Halton discrepancy...");
923
924 const Real * const discrepancy[8] = {
925 dim002DiscrRStartHalton, dim003DiscrRStartHalton,
926 dim005DiscrRStartHalton, dim010DiscrRStartHalton,
927 dim015DiscrRStartHalton, dim030DiscrRStartHalton,
928 dim050DiscrRStartHalton, dim100DiscrRStartHalton};
929
930 testGeneratorDiscrepancy(generatorFactory: HaltonFactory(true,false),
931 discrepancy,
932 "RandomStartHaltonDiscrepancy.txt",
933 "DiscrRStartHalton");
934}
935
936void LowDiscrepancyTest::testRandomShiftHaltonDiscrepancy() {
937
938 BOOST_TEST_MESSAGE("Testing random-shift Halton discrepancy...");
939
940 const Real * const discrepancy[8] = {
941 dim002DiscrRShiftHalton, dim003DiscrRShiftHalton,
942 dim005DiscrRShiftHalton, dim010DiscrRShiftHalton,
943 dim015DiscrRShiftHalton, dim030DiscrRShiftHalton,
944 dim050DiscrRShiftHalton, dim100DiscrRShiftHalton};
945
946 testGeneratorDiscrepancy(generatorFactory: HaltonFactory(false,true),
947 discrepancy,
948 "RandomShiftHaltonDiscrepancy.txt",
949 "DiscrRShiftHalton");
950}
951
952void LowDiscrepancyTest::testRandomStartRandomShiftHaltonDiscrepancy() {
953
954 BOOST_TEST_MESSAGE("Testing random-start, random-shift Halton discrepancy...");
955
956 const Real * const discrepancy[8] = {
957 dim002DiscrRStRShHalton, dim003DiscrRStRShHalton,
958 dim005DiscrRStRShHalton, dim010DiscrRStRShHalton,
959 dim015DiscrRStRShHalton, dim030DiscrRStRShHalton,
960 dim050DiscrRStRShHalton, dim100DiscrRStRShHalton};
961
962 testGeneratorDiscrepancy(generatorFactory: HaltonFactory(true,true),
963 discrepancy,
964 "RandomStartRandomShiftHaltonDiscrepancy.txt",
965 "DiscrRStRShHalton");
966}
967
968void LowDiscrepancyTest::testJackelSobolDiscrepancy() {
969
970 BOOST_TEST_MESSAGE("Testing Jaeckel-Sobol discrepancy...");
971
972 const Real * const discrepancy[8] = {
973 dim002Discr_Sobol, dim003Discr_Sobol,
974 dim005Discr_Sobol, dim010DiscrJackel_Sobol,
975 dim015DiscrJackel_Sobol, dim030DiscrJackel_Sobol,
976 dim050DiscrJackel_Sobol, dim100DiscrJackel_Sobol};
977
978 testGeneratorDiscrepancy(generatorFactory: SobolFactory(SobolRsg::Jaeckel),
979 discrepancy,
980 "JackelSobolDiscrepancy.txt",
981 "DiscrJackel_Sobol");
982}
983
984void LowDiscrepancyTest::testSobolLevitanSobolDiscrepancy() {
985
986 BOOST_TEST_MESSAGE("Testing Levitan-Sobol discrepancy...");
987
988 const Real * const discrepancy[8] = {
989 dim002Discr_Sobol, dim003Discr_Sobol,
990 dim005Discr_Sobol, dim010DiscrSobLev_Sobol,
991 dim015DiscrSobLev_Sobol, dim030DiscrSobLev_Sobol,
992 dim050DiscrSobLev_Sobol, dim100DiscrSobLev_Sobol};
993
994 testGeneratorDiscrepancy(generatorFactory: SobolFactory(SobolRsg::SobolLevitan),
995 discrepancy,
996 "SobolLevitanSobolDiscrepancy.txt",
997 "DiscrSobLev_Sobol");
998}
999
1000void LowDiscrepancyTest::testSobolLevitanLemieuxSobolDiscrepancy() {
1001
1002 BOOST_TEST_MESSAGE("Testing Levitan-Lemieux-Sobol discrepancy...");
1003
1004 const Real * const discrepancy[8] = {
1005 dim002Discr_Sobol, dim003Discr_Sobol,
1006 dim005Discr_Sobol, dim010DiscrSobLev_Sobol,
1007 dim015DiscrSobLev_Sobol, dim030DiscrSobLev_Sobol,
1008 dim050DiscrSobLem_Sobol, dim100DiscrSobLem_Sobol};
1009
1010 testGeneratorDiscrepancy(generatorFactory: SobolFactory(SobolRsg::SobolLevitanLemieux),
1011 discrepancy,
1012 "SobolLevitanLemieuxSobolDiscrepancy.txt",
1013 "DiscrSobLevLem_Sobol");
1014}
1015
1016void LowDiscrepancyTest::testUnitSobolDiscrepancy() {
1017
1018 BOOST_TEST_MESSAGE("Testing unit Sobol discrepancy...");
1019
1020 const Real* const discrepancy[8] = {dim002Discr_Unit_Sobol, dim003Discr_Unit_Sobol,
1021 dim005Discr_Unit_Sobol, dim010Discr_Unit_Sobol,
1022 dim015Discr_Unit_Sobol, dim030Discr_Unit_Sobol,
1023 dim050Discr_Unit_Sobol, dim100Discr_Unit_Sobol};
1024
1025 testGeneratorDiscrepancy(generatorFactory: SobolFactory(SobolRsg::Unit), discrepancy, "UnitSobolDiscrepancy.txt",
1026 "Discr__Unit_Sobol");
1027}
1028
1029
1030void LowDiscrepancyTest::testSobolSkipping() {
1031
1032 BOOST_TEST_MESSAGE("Testing Sobol sequence skipping...");
1033
1034 unsigned long seed = 42;
1035 Size dimensionality[] = { 1, 10, 100, 1000 };
1036 unsigned long skip[] = { 0, 1, 42, 512, 100000 };
1037 SobolRsg::DirectionIntegers integers[] = { SobolRsg::Unit,
1038 SobolRsg::Jaeckel,
1039 SobolRsg::SobolLevitan,
1040 SobolRsg::SobolLevitanLemieux };
1041
1042 for (auto& integer : integers) {
1043 for (Size& j : dimensionality) {
1044 for (unsigned long& k : skip) {
1045
1046 // extract n samples
1047 SobolRsg rsg1(j, seed, integer);
1048 for (Size l = 0; l < k; l++)
1049 rsg1.nextInt32Sequence();
1050
1051 // skip n samples at once
1052 SobolRsg rsg2(j, seed, integer);
1053 rsg2.skipTo(n: k);
1054
1055 // compare next 100 samples
1056 for (Size m = 0; m < 100; m++) {
1057 std::vector<std::uint_least32_t> s1 = rsg1.nextInt32Sequence();
1058 std::vector<std::uint_least32_t> s2 = rsg2.nextInt32Sequence();
1059 for (Size n = 0; n < s1.size(); n++) {
1060 if (s1[n] != s2[n]) {
1061 BOOST_ERROR("Mismatch after skipping:"
1062 << "\n size: " << j << "\n integers: " << integer
1063 << "\n skipped: " << k << "\n at index: " << n
1064 << "\n expected: " << s1[n] << "\n found: " << s2[n]);
1065 }
1066 }
1067 }
1068 }
1069 }
1070 }
1071}
1072
1073
1074test_suite* LowDiscrepancyTest::suite() {
1075 auto* suite = BOOST_TEST_SUITE("Low-discrepancy sequence tests");
1076
1077 suite->add(QUANTLIB_TEST_CASE(
1078 &LowDiscrepancyTest::testRandomizedLattices));
1079
1080
1081 suite->add(QUANTLIB_TEST_CASE(
1082 &LowDiscrepancyTest::testSeedGenerator));
1083
1084 suite->add(QUANTLIB_TEST_CASE(
1085 &LowDiscrepancyTest::testPolynomialsModuloTwo));
1086
1087 suite->add(QUANTLIB_TEST_CASE(&LowDiscrepancyTest::testSobol));
1088 suite->add(QUANTLIB_TEST_CASE(&LowDiscrepancyTest::testHalton));
1089 suite->add(QUANTLIB_TEST_CASE(&LowDiscrepancyTest::testFaure));
1090
1091 suite->add(QUANTLIB_TEST_CASE(
1092 &LowDiscrepancyTest::testMersenneTwisterDiscrepancy));
1093 suite->add(QUANTLIB_TEST_CASE(
1094 &LowDiscrepancyTest::testPlainHaltonDiscrepancy));
1095 suite->add(QUANTLIB_TEST_CASE(
1096 &LowDiscrepancyTest::testRandomStartHaltonDiscrepancy));
1097 suite->add(QUANTLIB_TEST_CASE(
1098 &LowDiscrepancyTest::testRandomShiftHaltonDiscrepancy));
1099 suite->add(QUANTLIB_TEST_CASE(
1100 &LowDiscrepancyTest::testRandomStartRandomShiftHaltonDiscrepancy));
1101 suite->add(QUANTLIB_TEST_CASE(
1102 &LowDiscrepancyTest::testUnitSobolDiscrepancy));
1103 suite->add(QUANTLIB_TEST_CASE(
1104 &LowDiscrepancyTest::testJackelSobolDiscrepancy));
1105 suite->add(QUANTLIB_TEST_CASE(
1106 &LowDiscrepancyTest::testSobolLevitanSobolDiscrepancy));
1107 suite->add(QUANTLIB_TEST_CASE(
1108 &LowDiscrepancyTest::testSobolLevitanLemieuxSobolDiscrepancy));
1109
1110 suite->add(QUANTLIB_TEST_CASE(&LowDiscrepancyTest::testSobolSkipping));
1111
1112 suite->add(QUANTLIB_TEST_CASE(
1113 &LowDiscrepancyTest::testRandomizedLowDiscrepancySequence));
1114
1115 return suite;
1116}
1117

source code of quantlib/test-suite/lowdiscrepancysequences.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.