QPR-12140 add unit test

Author: pcaspers

test-suite/lowdiscrepancysequences.cpp

@@ -30,6 +30,7 @@
 #include <ql/math/randomnumbers/randomizedlds.hpp>
 #include <ql/math/randomnumbers/randomsequencegenerator.hpp>
 #include <ql/math/randomnumbers/sobolrsg.hpp>
+#include <ql/math/randomnumbers/burley2020scrambledsobolrsg.hpp>
 #include <ql/utilities/dataformatters.hpp>
 #include <ql/math/randomnumbers/latticerules.hpp>
 #include <ql/math/randomnumbers/latticersg.hpp>
@@ -1070,6 +1071,62 @@ void LowDiscrepancyTest::testSobolSkipping() {
     }
 }
 
+void LowDiscrepancyTest::testHighDimensionalIntegrals() {
+
+    BOOST_TEST_MESSAGE("Testing High Dimensional Integrals");
+
+    /* We are running "Integration test 1, results for high dimensions" (Figure 9) from:
+
+       Sobol, Asotsky, Kreinin, Kucherenko: Construction and Comparison of High-Dimensional Sobol’
+       Generators, available at https://www.broda.co.uk/doc/HD_SobolGenerator.pdf
+
+       We check the error of Kuo1 (using Gray code and sequential numbers) roughly against what
+       their graph suggests. In addition we check the error of the Burley2020-scrambled version of
+       Kuo1 against what we experimentally find - the error turns out to be more than one order
+       better than the unscrambled version. */
+
+    auto integrand = [](const std::vector<Real>& c, const std::vector<Real>& x) {
+        Real p = 1.0;
+        for (Size i = 0; i < c.size(); ++i) {
+            p *= 1.0 + c[i] * (x[i] - 0.5);
+        }
+        return p;
+    };
+
+    Size N = 30031;
+
+    BOOST_TEST_MESSAGE("dimension,Sobol(Gray),Sobol(Seq),Burley2020");
+
+    std::vector<Size> dimension = {1000,2000,5000};
+    std::vector<std::vector<Real>> expectedOrderOfError = {
+        {-3.0, -3.0, -4.5}, {-2.5, -2.5, -4.0}, {-2.0, -2.0, -4.0}};
+
+    for (Size d = 0; d < dimension.size(); ++d) {
+
+        std::vector<Real> c1(dimension[d], 0.01);
+
+        SobolRsg s1(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, true);
+        SobolRsg s2(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, false);
+        Burley2020SobolRsg s3(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, 43);
+
+        Real I1 = 0.0, I2 = 0.0, I3 = 0.0;
+        for (Size i = 0; i < N; ++i) {
+            I1 += integrand(c1, s1.nextSequence().value) / static_cast<double>(N);
+            I2 += integrand(c1, s2.nextSequence().value) / static_cast<double>(N);
+            I3 += integrand(c1, s3.nextSequence().value) / static_cast<double>(N);
+        }
+
+        Real errOrder1 = std::log10(std::abs(I1-1.0));
+        Real errOrder2 = std::log10(std::abs(I2-1.0));
+        Real errOrder3 = std::log10(std::abs(I3-1.0));
+
+        BOOST_TEST_MESSAGE(dimension[d] << "," << errOrder1 << "," << errOrder2 << "," << errOrder3);
+
+        BOOST_CHECK(errOrder1 < expectedOrderOfError[d][0]);
+        BOOST_CHECK(errOrder2 < expectedOrderOfError[d][1]);
+        BOOST_CHECK(errOrder3 < expectedOrderOfError[d][2]);
+    }
+}
 
 test_suite* LowDiscrepancyTest::suite() {
     auto* suite = BOOST_TEST_SUITE("Low-discrepancy sequence tests");
@@ -1112,5 +1169,7 @@ test_suite* LowDiscrepancyTest::suite() {
     suite->add(QUANTLIB_TEST_CASE(
            &LowDiscrepancyTest::testRandomizedLowDiscrepancySequence));
 
+    suite->add(QUANTLIB_TEST_CASE(&LowDiscrepancyTest::testHighDimensionalIntegrals));
+
     return suite;
 }

test-suite/lowdiscrepancysequences.hpp

@@ -44,6 +44,7 @@ class LowDiscrepancyTest {
     static void testSobolLevitanSobolDiscrepancy();
     static void testSobolLevitanLemieuxSobolDiscrepancy();
     static void testUnitSobolDiscrepancy();
+    static void testHighDimensionalIntegrals();
 
     static void testRandomizedLowDiscrepancySequence();