libs/math/test/naive_monte_carlo_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2017
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7 #define BOOST_TEST_MODULE naive_monte_carlo_test
   8 #define BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
   9 #include <cmath>
  10 #include <ostream>
  11 #include <boost/lexical_cast.hpp>
  12 #include <boost/type_index.hpp>
  13 #include <boost/test/included/unit_test.hpp>
  14
  15 #include <boost/test/tools/floating_point_comparison.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17 #include <boost/math/quadrature/naive_monte_carlo.hpp>
  18
  19 using std::abs;
  20 using std::vector;
  21 using std::pair;
  22 using boost::math::constants::pi;
  23 using boost::math::quadrature::naive_monte_carlo;
  24
  25
  26 template<class Real>
  27 void test_pi_multithreaded()
  28 {
  29     std::cout << "Testing pi is calculated correctly (multithreaded) using Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  30     auto g = [](std::vector<Real> const & x)->Real {
  31         Real r = x[0]*x[0]+x[1]*x[1];
  32         if (r <= 1) {
  33           return 4;
  34         }
  35         return 0;
  36     };
  37
  38     std::vector<std::pair<Real, Real>> bounds{{Real(0), Real(1)}, {Real(0), Real(1)}};
  39     Real error_goal = 0.0002;
  40     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, error_goal,
  41                                           /*singular =*/ false,/* threads = */ 2, /* seed = */ 18012);
  42     auto task = mc.integrate();
  43     Real pi_estimated = task.get();
  44     if (abs(pi_estimated - pi<Real>())/pi<Real>() > 0.005) {
  45         std::cout << "Error in estimation of pi too high, function calls: " << mc.calls() << "\n";
  46         std::cout << "Final error estimate : " << mc.current_error_estimate() << "\n";
  47         std::cout << "Error goal           : " << error_goal << "\n";
  48         BOOST_CHECK_CLOSE_FRACTION(pi_estimated, pi<Real>(), 0.005);
  49     }
  50 }
  51
  52 template<class Real>
  53 void test_pi()
  54 {
  55     std::cout << "Testing pi is calculated correctly using Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  56     auto g = [](std::vector<Real> const & x)->Real
  57     {
  58         Real r = x[0]*x[0]+x[1]*x[1];
  59         if (r <= 1)
  60         {
  61             return 4;
  62         }
  63         return 0;
  64     };
  65
  66     std::vector<std::pair<Real, Real>> bounds{{Real(0), Real(1)}, {Real(0), Real(1)}};
  67     Real error_goal = (Real) 0.0002;
  68     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, error_goal,
  69                                             /*singular =*/ false,/* threads = */ 1, /* seed = */ 128402);
  70     auto task = mc.integrate();
  71     Real pi_estimated = task.get();
  72     if (abs(pi_estimated - pi<Real>())/pi<Real>() > 0.005)
  73     {
  74         std::cout << "Error in estimation of pi too high, function calls: " << mc.calls() << "\n";
  75         std::cout << "Final error estimate : " << mc.current_error_estimate() << "\n";
  76         std::cout << "Error goal           : " << error_goal << "\n";
  77         BOOST_CHECK_CLOSE_FRACTION(pi_estimated, pi<Real>(), 0.005);
  78     }
  79
  80 }
  81
  82 template<class Real>
  83 void test_constant()
  84 {
  85     std::cout << "Testing constants are integrated correctly using Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  86     auto g = [](std::vector<Real> const &)->Real
  87     {
  88       return 1;
  89     };
  90
  91     std::vector<std::pair<Real, Real>> bounds{{Real(0), Real(1)}, { Real(0), Real(1)}};
  92     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.0001,
  93                                             /* singular = */ false, /* threads = */ 1, /* seed = */ 87);
  94
  95     auto task = mc.integrate();
  96     Real one = task.get();
  97     BOOST_CHECK_CLOSE_FRACTION(one, 1, 0.001);
  98     BOOST_CHECK_SMALL(mc.current_error_estimate(), std::numeric_limits<Real>::epsilon());
  99     BOOST_CHECK(mc.calls() > 1000);
 100 }
 101
 102
 103 template<class Real>
 104 void test_exception_from_integrand()
 105 {
 106     std::cout << "Testing that a reasonable action is performed by the Monte-Carlo integrator when the integrand throws an exception on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 107     auto g = [](std::vector<Real> const & x)->Real
 108     {
 109         if (x[0] > 0.5 && x[0] < 0.5001)
 110         {
 111             throw std::domain_error("You have done something wrong.\n");
 112         }
 113         return (Real) 1;
 114     };
 115
 116     std::vector<std::pair<Real, Real>> bounds{{ Real(0), Real(1)}, { Real(0), Real(1)}};
 117     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.0001);
 118
 119     auto task = mc.integrate();
 120     bool caught_exception = false;
 121     try
 122     {
 123       Real result = task.get();
 124       // Get rid of unused variable warning:
 125       std::ostream cnull(0);
 126       cnull << result;
 127     }
 128     catch(std::exception const &)
 129     {
 130         caught_exception = true;
 131     }
 132     BOOST_CHECK(caught_exception);
 133 }
 134
 135
 136 template<class Real>
 137 void test_cancel_and_restart()
 138 {
 139     std::cout << "Testing that cancellation and restarting works on naive Monte-Carlo integration on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 140     Real exact = boost::lexical_cast<Real>("1.3932039296856768591842462603255");
 141     BOOST_CONSTEXPR const Real A = 1.0 / (pi<Real>() * pi<Real>() * pi<Real>());
 142     auto g = [&](std::vector<Real> const & x)->Real
 143     {
 144         return A / (1.0 - cos(x[0])*cos(x[1])*cos(x[2]));
 145     };
 146     vector<pair<Real, Real>> bounds{{ Real(0), pi<Real>()}, { Real(0), pi<Real>()}, { Real(0), pi<Real>()}};
 147     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.05, true, 1, 888889);
 148
 149     auto task = mc.integrate();
 150     mc.cancel();
 151     double y = task.get();
 152     // Super low tolerance; because it got canceled so fast:
 153     BOOST_CHECK_CLOSE_FRACTION(y, exact, 1.0);
 154
 155     mc.update_target_error((Real) 0.01);
 156     task = mc.integrate();
 157     y = task.get();
 158     BOOST_CHECK_CLOSE_FRACTION(y, exact, 0.1);
 159 }
 160
 161 template<class Real>
 162 void test_finite_singular_boundary()
 163 {
 164     std::cout << "Testing that finite singular boundaries work on naive Monte-Carlo integration on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 165     using std::pow;
 166     using std::log;
 167     auto g = [](std::vector<Real> const & x)->Real
 168     {
 169         // The first term is singular at x = 0.
 170         // The second at x = 1:
 171         return pow(log(1.0/x[0]), 2) + log1p(-x[0]);
 172     };
 173     vector<pair<Real, Real>> bounds{{Real(0), Real(1)}};
 174     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.01, true, 1, 1922);
 175
 176     auto task = mc.integrate();
 177
 178     double y = task.get();
 179     BOOST_CHECK_CLOSE_FRACTION(y, 1.0, 0.1);
 180 }
 181
 182 template<class Real>
 183 void test_multithreaded_variance()
 184 {
 185     std::cout << "Testing that variance computed by naive Monte-Carlo integration converges to integral formula on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 186     Real exact_variance = (Real) 1/(Real) 12;
 187     auto g = [&](std::vector<Real> const & x)->Real
 188     {
 189         return x[0];
 190     };
 191     vector<pair<Real, Real>> bounds{{ Real(0), Real(1)}};
 192     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, false, 2, 12341);
 193
 194     auto task = mc.integrate();
 195     Real y = task.get();
 196     BOOST_CHECK_CLOSE_FRACTION(y, 0.5, 0.01);
 197     BOOST_CHECK_CLOSE_FRACTION(mc.variance(), exact_variance, 0.05);
 198 }
 199
 200 template<class Real>
 201 void test_variance()
 202 {
 203     std::cout << "Testing that variance computed by naive Monte-Carlo integration converges to integral formula on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 204     Real exact_variance = (Real) 1/(Real) 12;
 205     auto g = [&](std::vector<Real> const & x)->Real
 206     {
 207         return x[0];
 208     };
 209     vector<pair<Real, Real>> bounds{{ Real(0), Real(1)}};
 210     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, false, 1, 12341);
 211
 212     auto task = mc.integrate();
 213     Real y = task.get();
 214     BOOST_CHECK_CLOSE_FRACTION(y, 0.5, 0.01);
 215     BOOST_CHECK_CLOSE_FRACTION(mc.variance(), exact_variance, 0.05);
 216 }
 217
 218 template<class Real, uint64_t dimension>
 219 void test_product()
 220 {
 221     std::cout << "Testing that product functions are integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 222     auto g = [&](std::vector<Real> const & x)->Real
 223     {
 224         double y = 1;
 225         for (uint64_t i = 0; i < x.size(); ++i)
 226         {
 227             y *= 2*x[i];
 228         }
 229         return y;
 230     };
 231
 232     vector<pair<Real, Real>> bounds(dimension);
 233     for (uint64_t i = 0; i < dimension; ++i)
 234     {
 235         bounds[i] = std::make_pair<Real, Real>(0, 1);
 236     }
 237     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, false, 1, 13999);
 238
 239     auto task = mc.integrate();
 240     Real y = task.get();
 241     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 242     using std::pow;
 243     Real exact_variance = pow(4.0/3.0, dimension) - 1;
 244     BOOST_CHECK_CLOSE_FRACTION(mc.variance(), exact_variance, 0.1);
 245 }
 246
 247 template<class Real, uint64_t dimension>
 248 void test_alternative_rng_1()
 249 {
 250     std::cout << "Testing that alternative RNGs work correctly using naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 251     auto g = [&](std::vector<Real> const & x)->Real
 252     {
 253         double y = 1;
 254         for (uint64_t i = 0; i < x.size(); ++i)
 255         {
 256             y *= 2*x[i];
 257         }
 258         return y;
 259     };
 260
 261     vector<pair<Real, Real>> bounds(dimension);
 262     for (uint64_t i = 0; i < dimension; ++i)
 263     {
 264         bounds[i] = std::make_pair<Real, Real>(0, 1);
 265     }
 266     std::cout << "Testing std::mt19937" << std::endl;
 267
 268     naive_monte_carlo<Real, decltype(g), std::mt19937> mc1(g, bounds, (Real) 0.001, false, 1, 1882);
 269
 270     auto task = mc1.integrate();
 271     Real y = task.get();
 272     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 273     using std::pow;
 274     Real exact_variance = pow(4.0/3.0, dimension) - 1;
 275     BOOST_CHECK_CLOSE_FRACTION(mc1.variance(), exact_variance, 0.05);
 276
 277     std::cout << "Testing std::knuth_b" << std::endl;
 278     naive_monte_carlo<Real, decltype(g), std::knuth_b> mc2(g, bounds, (Real) 0.001, false, 1, 1883);
 279     task = mc2.integrate();
 280     y = task.get();
 281     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 282
 283     std::cout << "Testing std::ranlux48" << std::endl;
 284     naive_monte_carlo<Real, decltype(g), std::ranlux48> mc3(g, bounds, (Real) 0.001, false, 1, 1884);
 285     task = mc3.integrate();
 286     y = task.get();
 287     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 288 }
 289
 290 template<class Real, uint64_t dimension>
 291 void test_alternative_rng_2()
 292 {
 293     std::cout << "Testing that alternative RNGs work correctly using naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 294     auto g = [&](std::vector<Real> const & x)->Real
 295     {
 296         double y = 1;
 297         for (uint64_t i = 0; i < x.size(); ++i)
 298         {
 299             y *= 2*x[i];
 300         }
 301         return y;
 302     };
 303
 304     vector<pair<Real, Real>> bounds(dimension);
 305     for (uint64_t i = 0; i < dimension; ++i)
 306     {
 307         bounds[i] = std::make_pair<Real, Real>(0, 1);
 308     }
 309
 310     std::cout << "Testing std::default_random_engine" << std::endl;
 311     naive_monte_carlo<Real, decltype(g), std::default_random_engine> mc4(g, bounds, (Real) 0.001, false, 1, 1884);
 312     auto task = mc4.integrate();
 313     Real y = task.get();
 314     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 315
 316     std::cout << "Testing std::minstd_rand" << std::endl;
 317     naive_monte_carlo<Real, decltype(g), std::minstd_rand> mc5(g, bounds, (Real) 0.001, false, 1, 1887);
 318     task = mc5.integrate();
 319     y = task.get();
 320     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 321
 322     std::cout << "Testing std::minstd_rand0" << std::endl;
 323     naive_monte_carlo<Real, decltype(g), std::minstd_rand0> mc6(g, bounds, (Real) 0.001, false, 1, 1889);
 324     task = mc6.integrate();
 325     y = task.get();
 326     BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 327
 328 }
 329
 330 template<class Real>
 331 void test_upper_bound_infinite()
 332 {
 333     std::cout << "Testing that infinite upper bounds are integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 334     auto g = [](std::vector<Real> const & x)->Real
 335     {
 336         return 1.0/(x[0]*x[0] + 1.0);
 337     };
 338
 339     vector<pair<Real, Real>> bounds(1);
 340     for (uint64_t i = 0; i < bounds.size(); ++i)
 341     {
 342         bounds[i] = std::make_pair<Real, Real>(0, std::numeric_limits<Real>::infinity());
 343     }
 344     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, true, 1, 8765);
 345     auto task = mc.integrate();
 346     Real y = task.get();
 347     BOOST_CHECK_CLOSE_FRACTION(y, boost::math::constants::half_pi<Real>(), 0.01);
 348 }
 349
 350 template<class Real>
 351 void test_lower_bound_infinite()
 352 {
 353     std::cout << "Testing that infinite lower bounds are integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 354     auto g = [](std::vector<Real> const & x)->Real
 355     {
 356         return 1.0/(x[0]*x[0] + 1.0);
 357     };
 358
 359     vector<pair<Real, Real>> bounds(1);
 360     for (uint64_t i = 0; i < bounds.size(); ++i)
 361     {
 362         bounds[i] = std::make_pair<Real, Real>(-std::numeric_limits<Real>::infinity(), 0);
 363     }
 364     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, true, 1, 1208);
 365
 366     auto task = mc.integrate();
 367     Real y = task.get();
 368     BOOST_CHECK_CLOSE_FRACTION(y, boost::math::constants::half_pi<Real>(), 0.01);
 369 }
 370
 371 template<class Real>
 372 void test_lower_bound_infinite2()
 373 {
 374     std::cout << "Testing that infinite lower bounds (2) are integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 375     auto g = [](std::vector<Real> const & x)->Real
 376     {
 377         // If x[0] = inf, this should blow up:
 378         return (x[0]*x[0])/(x[0]*x[0]*x[0]*x[0] + 1.0);
 379     };
 380
 381     vector<pair<Real, Real>> bounds(1);
 382     for (uint64_t i = 0; i < bounds.size(); ++i)
 383     {
 384         bounds[i] = std::make_pair<Real, Real>(-std::numeric_limits<Real>::infinity(), 0);
 385     }
 386     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, true, 1, 1208);
 387     auto task = mc.integrate();
 388     Real y = task.get();
 389     BOOST_CHECK_CLOSE_FRACTION(y, boost::math::constants::half_pi<Real>()/boost::math::constants::root_two<Real>(), 0.01);
 390 }
 391
 392 template<class Real>
 393 void test_double_infinite()
 394 {
 395     std::cout << "Testing that double infinite bounds are integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 396     auto g = [](std::vector<Real> const & x)->Real
 397     {
 398         return 1.0/(x[0]*x[0] + 1.0);
 399     };
 400
 401     vector<pair<Real, Real>> bounds(1);
 402     for (uint64_t i = 0; i < bounds.size(); ++i)
 403     {
 404         bounds[i] = std::make_pair<Real, Real>(-std::numeric_limits<Real>::infinity(), std::numeric_limits<Real>::infinity());
 405     }
 406     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, (Real) 0.001, true, 1, 1776);
 407
 408     auto task = mc.integrate();
 409     Real y = task.get();
 410     BOOST_CHECK_CLOSE_FRACTION(y, boost::math::constants::pi<Real>(), 0.01);
 411 }
 412
 413 template<class Real, uint64_t dimension>
 414 void test_radovic()
 415 {
 416     // See: Generalized Halton Sequences in 2008: A Comparative Study, function g1:
 417     std::cout << "Testing that the Radovic function is integrated correctly by naive Monte-Carlo on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 418     auto g = [](std::vector<Real> const & x)->Real
 419     {
 420         using std::abs;
 421         Real alpha = (Real)0.01;
 422         Real z = 1;
 423         for (uint64_t i = 0; i < dimension; ++i)
 424         {
 425             z *= (abs(4*x[i]-2) + alpha)/(1+alpha);
 426         }
 427         return z;
 428     };
 429
 430     vector<pair<Real, Real>> bounds(dimension);
 431     for (uint64_t i = 0; i < bounds.size(); ++i)
 432     {
 433         bounds[i] = std::make_pair<Real, Real>(0, 1);
 434     }
 435     Real error_goal = (Real) 0.001;
 436     naive_monte_carlo<Real, decltype(g)> mc(g, bounds, error_goal, false, 1, 1982);
 437
 438     auto task = mc.integrate();
 439     Real y = task.get();
 440     if (abs(y - 1) > 0.01)
 441     {
 442         std::cout << "Error in estimation of Radovic integral too high, function calls: " << mc.calls() << "\n";
 443         std::cout << "Final error estimate: " << mc.current_error_estimate() << std::endl;
 444         std::cout << "Error goal          : " << error_goal << std::endl;
 445         std::cout << "Variance estimate   : " << mc.variance() << std::endl;
 446         BOOST_CHECK_CLOSE_FRACTION(y, 1, 0.01);
 447     }
 448 }
 449
 450
 451 BOOST_AUTO_TEST_CASE(naive_monte_carlo_test)
 452 {
 453    std::cout << "Default hardware concurrency = " << std::thread::hardware_concurrency() << std::endl;
 454 #if !defined(TEST) || TEST == 1
 455     test_finite_singular_boundary<double>();
 456     test_finite_singular_boundary<float>();
 457 #endif
 458 #if !defined(TEST) || TEST == 2
 459     test_pi<float>();
 460     test_pi<double>();
 461 #endif
 462 #if !defined(TEST) || TEST == 3
 463     test_pi_multithreaded<float>();
 464     test_constant<float>();
 465 #endif
 466     //test_pi<long double>();
 467 #if !defined(TEST) || TEST == 4
 468     test_constant<double>();
 469     //test_constant<long double>();
 470     test_cancel_and_restart<float>();
 471 #endif
 472 #if !defined(TEST) || TEST == 5
 473     test_exception_from_integrand<float>();
 474     test_variance<float>();
 475 #endif
 476 #if !defined(TEST) || TEST == 6
 477     test_variance<double>();
 478     test_multithreaded_variance<double>();
 479 #endif
 480 #if !defined(TEST) || TEST == 7
 481     test_product<float, 1>();
 482     test_product<float, 2>();
 483 #endif
 484 #if !defined(TEST) || TEST == 8
 485     test_product<float, 3>();
 486     test_product<float, 4>();
 487     test_product<float, 5>();
 488 #endif
 489 #if !defined(TEST) || TEST == 9
 490     test_product<float, 6>();
 491     test_product<double, 1>();
 492 #endif
 493 #if !defined(TEST) || TEST == 10
 494     test_product<double, 2>();
 495 #endif
 496 #if !defined(TEST) || TEST == 11
 497     test_product<double, 3>();
 498     test_product<double, 4>();
 499 #endif
 500 #if !defined(TEST) || TEST == 12
 501     test_upper_bound_infinite<float>();
 502     test_upper_bound_infinite<double>();
 503 #endif
 504 #if !defined(TEST) || TEST == 13
 505     test_lower_bound_infinite<float>();
 506     test_lower_bound_infinite<double>();
 507 #endif
 508 #if !defined(TEST) || TEST == 14
 509     test_lower_bound_infinite2<float>();
 510 #endif
 511 #if !defined(TEST) || TEST == 15
 512     test_double_infinite<float>();
 513     test_double_infinite<double>();
 514 #endif
 515 #if !defined(TEST) || TEST == 16
 516     test_radovic<float, 1>();
 517     test_radovic<float, 2>();
 518 #endif
 519 #if !defined(TEST) || TEST == 17
 520     test_radovic<float, 3>();
 521     test_radovic<double, 1>();
 522 #endif
 523 #if !defined(TEST) || TEST == 18
 524     test_radovic<double, 2>();
 525     test_radovic<double, 3>();
 526 #endif
 527 #if !defined(TEST) || TEST == 19
 528     test_radovic<double, 4>();
 529     test_radovic<double, 5>();
 530 #endif
 531 #if !defined(TEST) || TEST == 20
 532     test_alternative_rng_1<float, 3>();
 533 #endif
 534 #if !defined(TEST) || TEST == 21
 535     test_alternative_rng_1<double, 3>();
 536 #endif
 537 #if !defined(TEST) || TEST == 22
 538     test_alternative_rng_2<float, 3>();
 539 #endif
 540 #if !defined(TEST) || TEST == 23
 541     test_alternative_rng_2<double, 3>();
 542 #endif
 543
 544 }