1[/ 2Copyright (c) 2018 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:naive_monte_carlo Naive Monte Carlo Integration] 9 10[heading Synopsis] 11 12 #include <boost/math/quadrature/naive_monte_carlo.hpp> 13 namespace boost { namespace math { namespace quadrature { 14 15 template<class Real, class F, class RNG = std::mt19937_64, class Policy = boost::math::policies::policy<>> 16 class naive_monte_carlo 17 { 18 public: 19 naive_monte_carlo(const F& integrand, 20 std::vector<std::pair<Real, Real>> const & bounds, 21 Real error_goal, 22 bool singular = true, 23 size_t threads = std::thread::hardware_concurrency()); 24 25 std::future<Real> integrate(); 26 27 void cancel(); 28 29 Real current_error_estimate() const; 30 31 std::chrono::duration<Real> estimated_time_to_completion() const; 32 33 void update_target_error(Real new_target_error); 34 35 Real progress() const; 36 37 Real current_estimate() const; 38 39 size_t calls() const; 40 }; 41 }}} // namespaces 42 43[heading Description] 44 45The class `naive_monte_carlo` performs Monte-Carlo integration on a square integrable function /f/ on a domain [Omega]. 46The theoretical background of Monte-Carlo integration is nicely discussed at [@https://en.wikipedia.org/wiki/Monte_Carlo_integration Wikipedia], 47and as such will not be discussed here. 48However, despite being "naive", 49it is a mistake to assume that naive Monte-Carlo integration is not powerful, 50as the simplicity of the method affords a robustness not easily provided by more sophisticated tools. 51The multithreaded nature of the routine allows us to compute a large number of sample points with great speed, 52and hence the slow convergence is mitigated by exploiting the full power of modern hardware. 53 54The naive Monte-Carlo integration provided by Boost exemplifies the programming techniques needed to cope with high-performance computing. 55For instance, since the convergence is only [bigo](N[super -1/2]), 56the compute time is very sensitive to the error goal. 57Users can easily specify an error goal which causes computation to last months-or just a few seconds. 58Without progress reporting, this situation is disorienting and causes the user to behave in a paranoid manner. 59Even with progress reporting, a user might need to cancel a job due to shifting priorities of the employing institution, 60and as such cancellation must be supported. 61A cancelled job which returns no results is wasted, so the cancellation must be graceful, 62returning the best estimate of the result thus far. 63In addition, a task might finish, and the user may well decide to try for a lower error bound. 64Hence restarting without loss of the preceding effort must be supported. 65Finally, on an HPC system, we generally wish to use all available threads. 66But if the computation is performed on a users workstation, 67employing every thread will cause the browser, email, or music apps to become unresponsive, 68so leaving a single thread available for other apps is appreciated. 69 70All these use cases are supported, so let's get to the code: 71 72 // Define a function to integrate: 73 auto g = [](std::vector<double> const & x) 74 { 75 constexpr const double A = 1.0 / (M_PI * M_PI * M_PI); 76 return A / (1.0 - cos(x[0])*cos(x[1])*cos(x[2])); 77 }; 78 std::vector<std::pair<double, double>> bounds{{0, M_PI}, {0, M_PI}, {0, M_PI}}; 79 double error_goal = 0.001 80 naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal); 81 82 std::future<double> task = mc.integrate(); 83 while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready) 84 { 85 // The user must decide on a reasonable way to display the progress depending on their environment: 86 display_progress(mc.progress(), 87 mc.current_error_estimate(), 88 mc.current_estimate(), 89 mc.estimated_time_to_completion()); 90 if (some_signal_heard()) 91 { 92 mc.cancel(); 93 std::cout << "\nCancelling because this is too slow!\n"; 94 } 95 } 96 double y = task.get(); 97 98First off, we define the function we wish to integrate. 99This function must accept a `std::vector<Real> const &`, and return a `Real`. 100Next, we define the domain of integration. 101Infinite domains are indicated by the bound `std::numeric_limits<Real>::infinity()`. 102The call 103 104 naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal); 105 106creates an instance of the monte carlo integrator. 107This is also where the number of threads can be set, for instance 108 109 naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal, true, std::thread::hardware_concurrency() - 1); 110 111might be more appropriate for running on a user's hardware (the default taking all the threads). 112The call to `integrate()` does not return the value of the integral, but rather a `std::future<Real>`. 113This allows us to do progress reporting from the master thread via 114 115 while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready) 116 { 117 // some reasonable method of displaying progress, based on the requirements of your app. 118 display_progress(mc.progress(), 119 mc.current_error_estimate(), 120 mc.current_estimate(), 121 mc.estimated_time_to_completion()); 122 } 123 124The file `example/naive_monte_carlo_example.cpp` has an implementation of `display_progress` which is reasonable for command line apps. 125In addition, we can call `mc.cancel()` in this loop to stop the integration. 126Progress reporting is especially useful if you accidentally pass in an integrand which is not square integrable-the variance increases without bound, and the progress decreases from some noisy initial value down to zero with time. 127Calling `task.get()` returns the current estimate. 128Once the future is ready, we can get the value of the integral via 129 130 double result = task.get(); 131 132At this point, the user may wish to reduce the error goal. 133This is achieved by 134 135 double new_target_error = 0.0005; 136 mc.update_target_error(new_target_error); 137 task = mc.integrate(); 138 y = task.get(); 139 140There is one additional "advanced" parameter: Whether or not the integrand is singular on the boundary. 141If the integrand is *not* singular on the boundary, then the integrand is evaluated over the closed set \u220F[sub i] [ /a/[sub /i/], /b/[sub /i/] ]. 142If the integrand is singular (the default) then the integrand is evaluated over the closed set \u220F[sub i] [ /a(1+[epsilon])/, /b(1-[epsilon])/ ]. 143(Note that there is sadly no such thing as an open set in floating point arithmetic.) 144When does the difference matter? 145Recall the stricture to never peel a high-dimensional orange, because when you do, nothing is left. 146The same idea applied here. 147The fraction of the volume within a distance [epsilon] of the boundary is approximately [epsilon]/d/, where /d/ is the number of dimensions. 148If the number of dimensions is large and the precision of the type is low, then it is possible that no correct digits will be obtained. 149If the integrand is singular on the boundary, you have no options; you simply must resort to higher precision computations. 150If the integrand is not singular on the boundary, then you can tell this to the integration routine via 151 152 naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal, /*singular = */ false); 153 154and this problem will not be encountered. 155In practice, you will need ~1,000 dimensions for this to be relevant in 16 bit floating point, ~100,000 dimensions in 32 bit floating point, 156and an astronomical number of dimensions in double precision. 157 158Finally, alternative random number generators may be provided to the class. 159The default random number generator is the standard library `std::mt19937_64`. 160However, here is an example which uses the 32-bit Mersenne twister random number generator instead: 161 162 naive_monte_carlo<Real, decltype(g), std::mt19937> mc(g, bounds, (Real) 0.001); 163 164[endsect] [/section:naive_monte_carlo Naive Monte Carlo Integration] 165 166