<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Naive Monte Carlo Integration</title>
<link rel="stylesheet" href="../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../index.html" title="Math Toolkit 3.0.0">
<link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation">
<link rel="prev" href="fourier_integrals.html" title="Fourier Integrals">
<link rel="next" href="wavelet_transforms.html" title="Wavelet Transforms">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
<td align="center"><a href="../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="fourier_integrals.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="wavelet_transforms.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="math_toolkit.naive_monte_carlo"></a><a class="link" href="naive_monte_carlo.html" title="Naive Monte Carlo Integration">Naive Monte Carlo Integration</a>
</h2></div></div></div>
<h4>
<a name="math_toolkit.naive_monte_carlo.h0"></a>
<span class="phrase"><a name="math_toolkit.naive_monte_carlo.synopsis"></a></span><a class="link" href="naive_monte_carlo.html#math_toolkit.naive_monte_carlo.synopsis">Synopsis</a>
</h4>
<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">quadrature</span><span class="special">/</span><span class="identifier">naive_monte_carlo</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
<span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">quadrature</span> <span class="special">{</span>
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">RNG</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">mt19937_64</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Policy</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special"><>></span>
<span class="keyword">class</span> <span class="identifier">naive_monte_carlo</span>
<span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
<span class="identifier">naive_monte_carlo</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span><span class="special">&</span> <span class="identifier">integrand</span><span class="special">,</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">>></span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">bounds</span><span class="special">,</span>
<span class="identifier">Real</span> <span class="identifier">error_goal</span><span class="special">,</span>
<span class="keyword">bool</span> <span class="identifier">singular</span> <span class="special">=</span> <span class="keyword">true</span><span class="special">,</span>
<span class="identifier">size_t</span> <span class="identifier">threads</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">thread</span><span class="special">::</span><span class="identifier">hardware_concurrency</span><span class="special">());</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">future</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="identifier">integrate</span><span class="special">();</span>
<span class="keyword">void</span> <span class="identifier">cancel</span><span class="special">();</span>
<span class="identifier">Real</span> <span class="identifier">current_error_estimate</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">chrono</span><span class="special">::</span><span class="identifier">duration</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="identifier">estimated_time_to_completion</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="keyword">void</span> <span class="identifier">update_target_error</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">new_target_error</span><span class="special">);</span>
<span class="identifier">Real</span> <span class="identifier">progress</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="identifier">Real</span> <span class="identifier">current_estimate</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="identifier">size_t</span> <span class="identifier">calls</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="special">};</span>
<span class="special">}}}</span> <span class="comment">// namespaces</span>
</pre>
<h4>
<a name="math_toolkit.naive_monte_carlo.h1"></a>
<span class="phrase"><a name="math_toolkit.naive_monte_carlo.description"></a></span><a class="link" href="naive_monte_carlo.html#math_toolkit.naive_monte_carlo.description">Description</a>
</h4>
<p>
The class <code class="computeroutput"><span class="identifier">naive_monte_carlo</span></code>
performs Monte-Carlo integration on a square integrable function <span class="emphasis"><em>f</em></span>
on a domain Ω. The theoretical background of Monte-Carlo integration is nicely
discussed at <a href="https://en.wikipedia.org/wiki/Monte_Carlo_integration" target="_top">Wikipedia</a>,
and as such will not be discussed here. However, despite being "naive",
it is a mistake to assume that naive Monte-Carlo integration is not powerful,
as the simplicity of the method affords a robustness not easily provided by
more sophisticated tools. The multithreaded nature of the routine allows us
to compute a large number of sample points with great speed, and hence the
slow convergence is mitigated by exploiting the full power of modern hardware.
</p>
<p>
The naive Monte-Carlo integration provided by Boost exemplifies the programming
techniques needed to cope with high-performance computing. For instance, since
the convergence is only 𝑶(N<sup>-1/2</sup>), the compute time is very sensitive to the
error goal. Users can easily specify an error goal which causes computation
to last months-or just a few seconds. Without progress reporting, this situation
is disorienting and causes the user to behave in a paranoid manner. Even with
progress reporting, a user might need to cancel a job due to shifting priorities
of the employing institution, and as such cancellation must be supported. A
cancelled job which returns no results is wasted, so the cancellation must
be graceful, returning the best estimate of the result thus far. In addition,
a task might finish, and the user may well decide to try for a lower error
bound. Hence restarting without loss of the preceding effort must be supported.
Finally, on an HPC system, we generally wish to use all available threads.
But if the computation is performed on a users workstation, employing every
thread will cause the browser, email, or music apps to become unresponsive,
so leaving a single thread available for other apps is appreciated.
</p>
<p>
All these use cases are supported, so let's get to the code:
</p>
<pre class="programlisting"><span class="comment">// Define a function to integrate:</span>
<span class="keyword">auto</span> <span class="identifier">g</span> <span class="special">=</span> <span class="special">[](</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">x</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">constexpr</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">A</span> <span class="special">=</span> <span class="number">1.0</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">M_PI</span> <span class="special">*</span> <span class="identifier">M_PI</span> <span class="special">*</span> <span class="identifier">M_PI</span><span class="special">);</span>
<span class="keyword">return</span> <span class="identifier">A</span> <span class="special">/</span> <span class="special">(</span><span class="number">1.0</span> <span class="special">-</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">])*</span><span class="identifier">cos</span><span class="special">(</span><span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">])*</span><span class="identifier">cos</span><span class="special">(</span><span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]));</span>
<span class="special">};</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">>></span> <span class="identifier">bounds</span><span class="special">{{</span><span class="number">0</span><span class="special">,</span> <span class="identifier">M_PI</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,</span> <span class="identifier">M_PI</span><span class="special">},</span> <span class="special">{</span><span class="number">0</span><span class="special">,</span> <span class="identifier">M_PI</span><span class="special">}};</span>
<span class="keyword">double</span> <span class="identifier">error_goal</span> <span class="special">=</span> <span class="number">0.001</span>
<span class="identifier">naive_monte_carlo</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(</span><span class="identifier">g</span><span class="special">)></span> <span class="identifier">mc</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">bounds</span><span class="special">,</span> <span class="identifier">error_goal</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">future</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">task</span> <span class="special">=</span> <span class="identifier">mc</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">();</span>
<span class="keyword">while</span> <span class="special">(</span><span class="identifier">task</span><span class="special">.</span><span class="identifier">wait_for</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">chrono</span><span class="special">::</span><span class="identifier">seconds</span><span class="special">(</span><span class="number">1</span><span class="special">))</span> <span class="special">!=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">future_status</span><span class="special">::</span><span class="identifier">ready</span><span class="special">)</span>
<span class="special">{</span>
<span class="comment">// The user must decide on a reasonable way to display the progress depending on their environment:</span>
<span class="identifier">display_progress</span><span class="special">(</span><span class="identifier">mc</span><span class="special">.</span><span class="identifier">progress</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_error_estimate</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_estimate</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">estimated_time_to_completion</span><span class="special">());</span>
<span class="keyword">if</span> <span class="special">(</span><span class="identifier">some_signal_heard</span><span class="special">())</span>
<span class="special">{</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">cancel</span><span class="special">();</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\nCancelling because this is too slow!\n"</span><span class="special">;</span>
<span class="special">}</span>
<span class="special">}</span>
<span class="keyword">double</span> <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">task</span><span class="special">.</span><span class="identifier">get</span><span class="special">();</span>
</pre>
<p>
First off, we define the function we wish to integrate. This function must
accept a <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="keyword">const</span> <span class="special">&</span></code>,
and return a <code class="computeroutput"><span class="identifier">Real</span></code>. Next, we
define the domain of integration. Infinite domains are indicated by the bound
<code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>::</span><span class="identifier">infinity</span><span class="special">()</span></code>.
The call
</p>
<pre class="programlisting"><span class="identifier">naive_monte_carlo</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(</span><span class="identifier">g</span><span class="special">)></span> <span class="identifier">mc</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">bounds</span><span class="special">,</span> <span class="identifier">error_goal</span><span class="special">);</span>
</pre>
<p>
creates an instance of the monte carlo integrator. This is also where the number
of threads can be set, for instance
</p>
<pre class="programlisting"><span class="identifier">naive_monte_carlo</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(</span><span class="identifier">g</span><span class="special">)></span> <span class="identifier">mc</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">bounds</span><span class="special">,</span> <span class="identifier">error_goal</span><span class="special">,</span> <span class="keyword">true</span><span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">thread</span><span class="special">::</span><span class="identifier">hardware_concurrency</span><span class="special">()</span> <span class="special">-</span> <span class="number">1</span><span class="special">);</span>
</pre>
<p>
might be more appropriate for running on a user's hardware (the default taking
all the threads). The call to <code class="computeroutput"><span class="identifier">integrate</span><span class="special">()</span></code> does not return the value of the integral,
but rather a <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">future</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span></code>.
This allows us to do progress reporting from the master thread via
</p>
<pre class="programlisting"><span class="keyword">while</span> <span class="special">(</span><span class="identifier">task</span><span class="special">.</span><span class="identifier">wait_for</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">chrono</span><span class="special">::</span><span class="identifier">seconds</span><span class="special">(</span><span class="number">1</span><span class="special">))</span> <span class="special">!=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">future_status</span><span class="special">::</span><span class="identifier">ready</span><span class="special">)</span>
<span class="special">{</span>
<span class="comment">// some reasonable method of displaying progress, based on the requirements of your app.</span>
<span class="identifier">display_progress</span><span class="special">(</span><span class="identifier">mc</span><span class="special">.</span><span class="identifier">progress</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_error_estimate</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_estimate</span><span class="special">(),</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">estimated_time_to_completion</span><span class="special">());</span>
<span class="special">}</span>
</pre>
<p>
The file <code class="computeroutput"><span class="identifier">example</span><span class="special">/</span><span class="identifier">naive_monte_carlo_example</span><span class="special">.</span><span class="identifier">cpp</span></code> has an implementation of <code class="computeroutput"><span class="identifier">display_progress</span></code> which is reasonable for
command line apps. In addition, we can call <code class="computeroutput"><span class="identifier">mc</span><span class="special">.</span><span class="identifier">cancel</span><span class="special">()</span></code>
in this loop to stop the integration. 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. Calling <code class="computeroutput"><span class="identifier">task</span><span class="special">.</span><span class="identifier">get</span><span class="special">()</span></code>
returns the current estimate. Once the future is ready, we can get the value
of the integral via
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">task</span><span class="special">.</span><span class="identifier">get</span><span class="special">();</span>
</pre>
<p>
At this point, the user may wish to reduce the error goal. This is achieved
by
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">new_target_error</span> <span class="special">=</span> <span class="number">0.0005</span><span class="special">;</span>
<span class="identifier">mc</span><span class="special">.</span><span class="identifier">update_target_error</span><span class="special">(</span><span class="identifier">new_target_error</span><span class="special">);</span>
<span class="identifier">task</span> <span class="special">=</span> <span class="identifier">mc</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">();</span>
<span class="identifier">y</span> <span class="special">=</span> <span class="identifier">task</span><span class="special">.</span><span class="identifier">get</span><span class="special">();</span>
</pre>
<p>
There is one additional "advanced" parameter: Whether or not the
integrand is singular on the boundary. If the integrand is <span class="bold"><strong>not</strong></span>
singular on the boundary, then the integrand is evaluated over the closed set
∏<sub>i</sub> [ <span class="emphasis"><em>a</em></span><sub><span class="emphasis"><em>i</em></span></sub>, <span class="emphasis"><em>b</em></span><sub><span class="emphasis"><em>i</em></span></sub> ].
If the integrand is singular (the default) then the integrand is evaluated
over the closed set ∏<sub>i</sub> [ /a(1+ε)/, /b(1-ε)/ ]. (Note that there is sadly
no such thing as an open set in floating point arithmetic.) When does the difference
matter? Recall the stricture to never peel a high-dimensional orange, because
when you do, nothing is left. The same idea applied here. The fraction of the
volume within a distance ε of the boundary is approximately ε<span class="emphasis"><em>d</em></span>,
where <span class="emphasis"><em>d</em></span> is the number of dimensions. 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. If the integrand is singular on the
boundary, you have no options; you simply must resort to higher precision computations.
If the integrand is not singular on the boundary, then you can tell this to
the integration routine via
</p>
<pre class="programlisting"><span class="identifier">naive_monte_carlo</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(</span><span class="identifier">g</span><span class="special">)></span> <span class="identifier">mc</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">bounds</span><span class="special">,</span> <span class="identifier">error_goal</span><span class="special">,</span> <span class="comment">/*singular = */</span> <span class="keyword">false</span><span class="special">);</span>
</pre>
<p>
and this problem will not be encountered. 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, and an astronomical number of dimensions in double
precision.
</p>
<p>
Finally, alternative random number generators may be provided to the class.
The default random number generator is the standard library <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">mt19937_64</span></code>.
However, here is an example which uses the 32-bit Mersenne twister random number
generator instead:
</p>
<pre class="programlisting"><span class="identifier">naive_monte_carlo</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="keyword">decltype</span><span class="special">(</span><span class="identifier">g</span><span class="special">),</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">mt19937</span><span class="special">></span> <span class="identifier">mc</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">bounds</span><span class="special">,</span> <span class="special">(</span><span class="identifier">Real</span><span class="special">)</span> <span class="number">0.001</span><span class="special">);</span>
</pre>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="fourier_integrals.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="wavelet_transforms.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>