libs/math/doc/quadrature/naive_monte_carlo.qbk

   1 [/
   2 Copyright (c) 2018 Nick Thompson
   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
   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
  45 The class `naive_monte_carlo` performs Monte-Carlo integration on a square integrable function /f/ on a domain [Omega].
  46 The theoretical background of Monte-Carlo integration is nicely discussed at [@https://en.wikipedia.org/wiki/Monte_Carlo_integration Wikipedia],
  47 and as such will not be discussed here.
  48 However, despite being "naive",
  49 it is a mistake to assume that naive Monte-Carlo integration is not powerful,
  50 as the simplicity of the method affords a robustness not easily provided by more sophisticated tools.
  51 The multithreaded nature of the routine allows us to compute a large number of sample points with great speed,
  52 and hence the slow convergence is mitigated by exploiting the full power of modern hardware.
  53
  54 The naive Monte-Carlo integration provided by Boost exemplifies the programming techniques needed to cope with high-performance computing.
  55 For instance, since the convergence is only [bigo](N[super -1/2]),
  56 the compute time is very sensitive to the error goal.
  57 Users can easily specify an error goal which causes computation to last months-or just a few seconds.
  58 Without progress reporting, this situation is disorienting and causes the user to behave in a paranoid manner.
  59 Even with progress reporting, a user might need to cancel a job due to shifting priorities of the employing institution,
  60 and as such cancellation must be supported.
  61 A cancelled job which returns no results is wasted, so the cancellation must be graceful,
  62 returning the best estimate of the result thus far.
  63 In addition, a task might finish, and the user may well decide to try for a lower error bound.
  64 Hence restarting without loss of the preceding effort must be supported.
  65 Finally, on an HPC system, we generally wish to use all available threads.
  66 But if the computation is performed on a users workstation,
  67 employing every thread will cause the browser, email, or music apps to become unresponsive,
  68 so leaving a single thread available for other apps is appreciated.
  69
  70 All 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
  98 First off, we define the function we wish to integrate.
  99 This function must accept a `std::vector<Real> const &`, and return a `Real`.
 100 Next, we define the domain of integration.
 101 Infinite domains are indicated by the bound `std::numeric_limits<Real>::infinity()`.
 102 The call
 103
 104     naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal);
 105
 106 creates an instance of the monte carlo integrator.
 107 This 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
 111 might be more appropriate for running on a user's hardware (the default taking all the threads).
 112 The call to `integrate()` does not return the value of the integral, but rather a `std::future<Real>`.
 113 This 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
 124 The file `example/naive_monte_carlo_example.cpp` has an implementation of `display_progress` which is reasonable for command line apps.
 125 In addition, we can call `mc.cancel()` in this loop to stop the integration.
 126 Progress 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.
 127 Calling `task.get()` returns the current estimate.
 128 Once the future is ready, we can get the value of the integral via
 129
 130     double result = task.get();
 131
 132 At this point, the user may wish to reduce the error goal.
 133 This 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
 140 There is one additional "advanced" parameter: Whether or not the integrand is singular on the boundary.
 141 If 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/] ].
 142 If 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.)
 144 When does the difference matter?
 145 Recall the stricture to never peel a high-dimensional orange, because when you do, nothing is left.
 146 The same idea applied here.
 147 The fraction of the volume within a distance [epsilon] of the boundary is approximately [epsilon]/d/, where /d/ is the number of dimensions.
 148 If 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.
 149 If the integrand is singular on the boundary, you have no options; you simply must resort to higher precision computations.
 150 If 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
 154 and this problem will not be encountered.
 155 In 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,
 156 and an astronomical number of dimensions in double precision.
 157
 158 Finally, alternative random number generators may be provided to the class.
 159 The default random number generator is the standard library `std::mt19937_64`.
 160 However, 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