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