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)
7 #define BOOST_TEST_MODULE naive_monte_carlo_test
8 #define BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
11 #include <boost/lexical_cast.hpp>
12 #include <boost/type_index.hpp>
13 #include <boost/test/included/unit_test.hpp>
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>
22 using boost::math::constants::pi
;
23 using boost::math::quadrature::naive_monte_carlo
;
27 void test_pi_multithreaded()
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];
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);
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
58 Real r
= x
[0]*x
[0]+x
[1]*x
[1];
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)
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);
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
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);
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);
104 void test_exception_from_integrand()
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
109 if (x
[0] > 0.5 && x
[0] < 0.5001)
111 throw std::domain_error("You have done something wrong.\n");
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);
119 auto task
= mc
.integrate();
120 bool caught_exception
= false;
123 Real result
= task
.get();
124 // Get rid of unused variable warning:
125 std::ostream
cnull(0);
128 catch(std::exception
const &)
130 caught_exception
= true;
132 BOOST_CHECK(caught_exception
);
137 void test_cancel_and_restart()
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
144 return A
/ (1.0 - cos(x
[0])*cos(x
[1])*cos(x
[2]));
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);
149 auto task
= mc
.integrate();
151 double y
= task
.get();
152 // Super low tolerance; because it got canceled so fast:
153 BOOST_CHECK_CLOSE_FRACTION(y
, exact
, 1.0);
155 mc
.update_target_error((Real
) 0.01);
156 task
= mc
.integrate();
158 BOOST_CHECK_CLOSE_FRACTION(y
, exact
, 0.1);
162 void test_finite_singular_boundary()
164 std::cout
<< "Testing that finite singular boundaries work on naive Monte-Carlo integration on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
167 auto g
= [](std::vector
<Real
> const & x
)->Real
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]);
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);
176 auto task
= mc
.integrate();
178 double y
= task
.get();
179 BOOST_CHECK_CLOSE_FRACTION(y
, 1.0, 0.1);
183 void test_multithreaded_variance()
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
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);
194 auto task
= mc
.integrate();
196 BOOST_CHECK_CLOSE_FRACTION(y
, 0.5, 0.01);
197 BOOST_CHECK_CLOSE_FRACTION(mc
.variance(), exact_variance
, 0.05);
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
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);
212 auto task
= mc
.integrate();
214 BOOST_CHECK_CLOSE_FRACTION(y
, 0.5, 0.01);
215 BOOST_CHECK_CLOSE_FRACTION(mc
.variance(), exact_variance
, 0.05);
218 template<class Real
, uint64_t dimension
>
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
225 for (uint64_t i
= 0; i
< x
.size(); ++i
)
232 vector
<pair
<Real
, Real
>> bounds(dimension
);
233 for (uint64_t i
= 0; i
< dimension
; ++i
)
235 bounds
[i
] = std::make_pair
<Real
, Real
>(0, 1);
237 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, (Real
) 0.001, false, 1, 13999);
239 auto task
= mc
.integrate();
241 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
243 Real exact_variance
= pow(4.0/3.0, dimension
) - 1;
244 BOOST_CHECK_CLOSE_FRACTION(mc
.variance(), exact_variance
, 0.1);
247 template<class Real
, uint64_t dimension
>
248 void test_alternative_rng_1()
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
254 for (uint64_t i
= 0; i
< x
.size(); ++i
)
261 vector
<pair
<Real
, Real
>> bounds(dimension
);
262 for (uint64_t i
= 0; i
< dimension
; ++i
)
264 bounds
[i
] = std::make_pair
<Real
, Real
>(0, 1);
266 std::cout
<< "Testing std::mt19937" << std::endl
;
268 naive_monte_carlo
<Real
, decltype(g
), std::mt19937
> mc1(g
, bounds
, (Real
) 0.001, false, 1, 1882);
270 auto task
= mc1
.integrate();
272 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
274 Real exact_variance
= pow(4.0/3.0, dimension
) - 1;
275 BOOST_CHECK_CLOSE_FRACTION(mc1
.variance(), exact_variance
, 0.05);
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();
281 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
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();
287 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
290 template<class Real
, uint64_t dimension
>
291 void test_alternative_rng_2()
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
297 for (uint64_t i
= 0; i
< x
.size(); ++i
)
304 vector
<pair
<Real
, Real
>> bounds(dimension
);
305 for (uint64_t i
= 0; i
< dimension
; ++i
)
307 bounds
[i
] = std::make_pair
<Real
, Real
>(0, 1);
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();
314 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
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();
320 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
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();
326 BOOST_CHECK_CLOSE_FRACTION(y
, 1, 0.01);
331 void test_upper_bound_infinite()
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
336 return 1.0/(x
[0]*x
[0] + 1.0);
339 vector
<pair
<Real
, Real
>> bounds(1);
340 for (uint64_t i
= 0; i
< bounds
.size(); ++i
)
342 bounds
[i
] = std::make_pair
<Real
, Real
>(0, std::numeric_limits
<Real
>::infinity());
344 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, (Real
) 0.001, true, 1, 8765);
345 auto task
= mc
.integrate();
347 BOOST_CHECK_CLOSE_FRACTION(y
, boost::math::constants::half_pi
<Real
>(), 0.01);
351 void test_lower_bound_infinite()
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
356 return 1.0/(x
[0]*x
[0] + 1.0);
359 vector
<pair
<Real
, Real
>> bounds(1);
360 for (uint64_t i
= 0; i
< bounds
.size(); ++i
)
362 bounds
[i
] = std::make_pair
<Real
, Real
>(-std::numeric_limits
<Real
>::infinity(), 0);
364 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, (Real
) 0.001, true, 1, 1208);
366 auto task
= mc
.integrate();
368 BOOST_CHECK_CLOSE_FRACTION(y
, boost::math::constants::half_pi
<Real
>(), 0.01);
372 void test_lower_bound_infinite2()
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
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);
381 vector
<pair
<Real
, Real
>> bounds(1);
382 for (uint64_t i
= 0; i
< bounds
.size(); ++i
)
384 bounds
[i
] = std::make_pair
<Real
, Real
>(-std::numeric_limits
<Real
>::infinity(), 0);
386 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, (Real
) 0.001, true, 1, 1208);
387 auto task
= mc
.integrate();
389 BOOST_CHECK_CLOSE_FRACTION(y
, boost::math::constants::half_pi
<Real
>()/boost::math::constants::root_two
<Real
>(), 0.01);
393 void test_double_infinite()
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
398 return 1.0/(x
[0]*x
[0] + 1.0);
401 vector
<pair
<Real
, Real
>> bounds(1);
402 for (uint64_t i
= 0; i
< bounds
.size(); ++i
)
404 bounds
[i
] = std::make_pair
<Real
, Real
>(-std::numeric_limits
<Real
>::infinity(), std::numeric_limits
<Real
>::infinity());
406 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, (Real
) 0.001, true, 1, 1776);
408 auto task
= mc
.integrate();
410 BOOST_CHECK_CLOSE_FRACTION(y
, boost::math::constants::pi
<Real
>(), 0.01);
413 template<class Real
, uint64_t dimension
>
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
421 Real alpha
= (Real
)0.01;
423 for (uint64_t i
= 0; i
< dimension
; ++i
)
425 z
*= (abs(4*x
[i
]-2) + alpha
)/(1+alpha
);
430 vector
<pair
<Real
, Real
>> bounds(dimension
);
431 for (uint64_t i
= 0; i
< bounds
.size(); ++i
)
433 bounds
[i
] = std::make_pair
<Real
, Real
>(0, 1);
435 Real error_goal
= (Real
) 0.001;
436 naive_monte_carlo
<Real
, decltype(g
)> mc(g
, bounds
, error_goal
, false, 1, 1982);
438 auto task
= mc
.integrate();
440 if (abs(y
- 1) > 0.01)
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);
451 BOOST_AUTO_TEST_CASE(naive_monte_carlo_test
)
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>();
458 #if !defined(TEST) || TEST == 2
462 #if !defined(TEST) || TEST == 3
463 test_pi_multithreaded
<float>();
464 test_constant
<float>();
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>();
472 #if !defined(TEST) || TEST == 5
473 test_exception_from_integrand
<float>();
474 test_variance
<float>();
476 #if !defined(TEST) || TEST == 6
477 test_variance
<double>();
478 test_multithreaded_variance
<double>();
480 #if !defined(TEST) || TEST == 7
481 test_product
<float, 1>();
482 test_product
<float, 2>();
484 #if !defined(TEST) || TEST == 8
485 test_product
<float, 3>();
486 test_product
<float, 4>();
487 test_product
<float, 5>();
489 #if !defined(TEST) || TEST == 9
490 test_product
<float, 6>();
491 test_product
<double, 1>();
493 #if !defined(TEST) || TEST == 10
494 test_product
<double, 2>();
496 #if !defined(TEST) || TEST == 11
497 test_product
<double, 3>();
498 test_product
<double, 4>();
500 #if !defined(TEST) || TEST == 12
501 test_upper_bound_infinite
<float>();
502 test_upper_bound_infinite
<double>();
504 #if !defined(TEST) || TEST == 13
505 test_lower_bound_infinite
<float>();
506 test_lower_bound_infinite
<double>();
508 #if !defined(TEST) || TEST == 14
509 test_lower_bound_infinite2
<float>();
511 #if !defined(TEST) || TEST == 15
512 test_double_infinite
<float>();
513 test_double_infinite
<double>();
515 #if !defined(TEST) || TEST == 16
516 test_radovic
<float, 1>();
517 test_radovic
<float, 2>();
519 #if !defined(TEST) || TEST == 17
520 test_radovic
<float, 3>();
521 test_radovic
<double, 1>();
523 #if !defined(TEST) || TEST == 18
524 test_radovic
<double, 2>();
525 test_radovic
<double, 3>();
527 #if !defined(TEST) || TEST == 19
528 test_radovic
<double, 4>();
529 test_radovic
<double, 5>();
531 #if !defined(TEST) || TEST == 20
532 test_alternative_rng_1
<float, 3>();
534 #if !defined(TEST) || TEST == 21
535 test_alternative_rng_1
<double, 3>();
537 #if !defined(TEST) || TEST == 22
538 test_alternative_rng_2
<float, 3>();
540 #if !defined(TEST) || TEST == 23
541 test_alternative_rng_2
<double, 3>();