1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Naive Monte Carlo Integration</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation"> 9<link rel="prev" href="fourier_integrals.html" title="Fourier Integrals"> 10<link rel="next" href="wavelet_transforms.html" title="Wavelet Transforms"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<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> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<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> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.naive_monte_carlo.h0"></a> 31 <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> 32 </h4> 33<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> 34<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> 35 36<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> 37<span class="keyword">class</span> <span class="identifier">naive_monte_carlo</span> 38<span class="special">{</span> 39<span class="keyword">public</span><span class="special">:</span> 40 <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> 41 <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> 42 <span class="identifier">Real</span> <span class="identifier">error_goal</span><span class="special">,</span> 43 <span class="keyword">bool</span> <span class="identifier">singular</span> <span class="special">=</span> <span class="keyword">true</span><span class="special">,</span> 44 <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> 45 46 <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> 47 48 <span class="keyword">void</span> <span class="identifier">cancel</span><span class="special">();</span> 49 50 <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> 51 52 <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> 53 54 <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> 55 56 <span class="identifier">Real</span> <span class="identifier">progress</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 57 58 <span class="identifier">Real</span> <span class="identifier">current_estimate</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 59 60 <span class="identifier">size_t</span> <span class="identifier">calls</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 61<span class="special">};</span> 62<span class="special">}}}</span> <span class="comment">// namespaces</span> 63</pre> 64<h4> 65<a name="math_toolkit.naive_monte_carlo.h1"></a> 66 <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> 67 </h4> 68<p> 69 The class <code class="computeroutput"><span class="identifier">naive_monte_carlo</span></code> 70 performs Monte-Carlo integration on a square integrable function <span class="emphasis"><em>f</em></span> 71 on a domain Ω. The theoretical background of Monte-Carlo integration is nicely 72 discussed at <a href="https://en.wikipedia.org/wiki/Monte_Carlo_integration" target="_top">Wikipedia</a>, 73 and as such will not be discussed here. However, despite being "naive", 74 it is a mistake to assume that naive Monte-Carlo integration is not powerful, 75 as the simplicity of the method affords a robustness not easily provided by 76 more sophisticated tools. The multithreaded nature of the routine allows us 77 to compute a large number of sample points with great speed, and hence the 78 slow convergence is mitigated by exploiting the full power of modern hardware. 79 </p> 80<p> 81 The naive Monte-Carlo integration provided by Boost exemplifies the programming 82 techniques needed to cope with high-performance computing. For instance, since 83 the convergence is only ��(N<sup>-1/2</sup>), the compute time is very sensitive to the 84 error goal. Users can easily specify an error goal which causes computation 85 to last months-or just a few seconds. Without progress reporting, this situation 86 is disorienting and causes the user to behave in a paranoid manner. Even with 87 progress reporting, a user might need to cancel a job due to shifting priorities 88 of the employing institution, and as such cancellation must be supported. A 89 cancelled job which returns no results is wasted, so the cancellation must 90 be graceful, returning the best estimate of the result thus far. In addition, 91 a task might finish, and the user may well decide to try for a lower error 92 bound. Hence restarting without loss of the preceding effort must be supported. 93 Finally, on an HPC system, we generally wish to use all available threads. 94 But if the computation is performed on a users workstation, employing every 95 thread will cause the browser, email, or music apps to become unresponsive, 96 so leaving a single thread available for other apps is appreciated. 97 </p> 98<p> 99 All these use cases are supported, so let's get to the code: 100 </p> 101<pre class="programlisting"><span class="comment">// Define a function to integrate:</span> 102<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> 103<span class="special">{</span> 104 <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> 105 <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> 106<span class="special">};</span> 107<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> 108<span class="keyword">double</span> <span class="identifier">error_goal</span> <span class="special">=</span> <span class="number">0.001</span> 109<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> 110 111<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> 112<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> 113<span class="special">{</span> 114 <span class="comment">// The user must decide on a reasonable way to display the progress depending on their environment:</span> 115 <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> 116 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_error_estimate</span><span class="special">(),</span> 117 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_estimate</span><span class="special">(),</span> 118 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">estimated_time_to_completion</span><span class="special">());</span> 119 <span class="keyword">if</span> <span class="special">(</span><span class="identifier">some_signal_heard</span><span class="special">())</span> 120 <span class="special">{</span> 121 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">cancel</span><span class="special">();</span> 122 <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> 123 <span class="special">}</span> 124<span class="special">}</span> 125<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> 126</pre> 127<p> 128 First off, we define the function we wish to integrate. This function must 129 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>, 130 and return a <code class="computeroutput"><span class="identifier">Real</span></code>. Next, we 131 define the domain of integration. Infinite domains are indicated by the bound 132 <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>. 133 The call 134 </p> 135<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> 136</pre> 137<p> 138 creates an instance of the monte carlo integrator. This is also where the number 139 of threads can be set, for instance 140 </p> 141<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> 142</pre> 143<p> 144 might be more appropriate for running on a user's hardware (the default taking 145 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, 146 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>. 147 This allows us to do progress reporting from the master thread via 148 </p> 149<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> 150<span class="special">{</span> 151 <span class="comment">// some reasonable method of displaying progress, based on the requirements of your app.</span> 152 <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> 153 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_error_estimate</span><span class="special">(),</span> 154 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">current_estimate</span><span class="special">(),</span> 155 <span class="identifier">mc</span><span class="special">.</span><span class="identifier">estimated_time_to_completion</span><span class="special">());</span> 156<span class="special">}</span> 157</pre> 158<p> 159 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 160 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> 161 in this loop to stop the integration. Progress reporting is especially useful 162 if you accidentally pass in an integrand which is not square integrable-the 163 variance increases without bound, and the progress decreases from some noisy 164 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> 165 returns the current estimate. Once the future is ready, we can get the value 166 of the integral via 167 </p> 168<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> 169</pre> 170<p> 171 At this point, the user may wish to reduce the error goal. This is achieved 172 by 173 </p> 174<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> 175<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> 176<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> 177<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> 178</pre> 179<p> 180 There is one additional "advanced" parameter: Whether or not the 181 integrand is singular on the boundary. If the integrand is <span class="bold"><strong>not</strong></span> 182 singular on the boundary, then the integrand is evaluated over the closed set 183 ∏<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> ]. 184 If the integrand is singular (the default) then the integrand is evaluated 185 over the closed set ∏<sub>i</sub> [ /a(1+ε)/, /b(1-ε)/ ]. (Note that there is sadly 186 no such thing as an open set in floating point arithmetic.) When does the difference 187 matter? Recall the stricture to never peel a high-dimensional orange, because 188 when you do, nothing is left. The same idea applied here. The fraction of the 189 volume within a distance ε of the boundary is approximately ε<span class="emphasis"><em>d</em></span>, 190 where <span class="emphasis"><em>d</em></span> is the number of dimensions. If the number of 191 dimensions is large and the precision of the type is low, then it is possible 192 that no correct digits will be obtained. If the integrand is singular on the 193 boundary, you have no options; you simply must resort to higher precision computations. 194 If the integrand is not singular on the boundary, then you can tell this to 195 the integration routine via 196 </p> 197<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> 198</pre> 199<p> 200 and this problem will not be encountered. In practice, you will need ~1,000 201 dimensions for this to be relevant in 16 bit floating point, ~100,000 dimensions 202 in 32 bit floating point, and an astronomical number of dimensions in double 203 precision. 204 </p> 205<p> 206 Finally, alternative random number generators may be provided to the class. 207 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>. 208 However, here is an example which uses the 32-bit Mersenne twister random number 209 generator instead: 210 </p> 211<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> 212</pre> 213</div> 214<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 215<td align="left"></td> 216<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 217 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 218 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 219 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 220 Daryle Walker and Xiaogang Zhang<p> 221 Distributed under the Boost Software License, Version 1.0. (See accompanying 222 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>) 223 </p> 224</div></td> 225</tr></table> 226<hr> 227<div class="spirit-nav"> 228<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> 229</div> 230</body> 231</html> 232
