devn00b
/
EQ2EMu


			
				
					
						
						
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231
							<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<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 2.11.0">
<link rel="up" href="../quadrature.html" title="Chapter&#160;13.&#160;Quadrature and Differentiation">
<link rel="prev" href="fourier_integrals.html" title="Fourier Integrals">
<link rel="next" href="diff.html" title="Numerical Differentiation">
</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="diff.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">&lt;</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">&gt;</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">&lt;</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">&lt;&gt;&gt;</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">&amp;</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">&lt;</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">&gt;&gt;</span> <span class="keyword">const</span> <span class="special">&amp;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;</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 &#937;. 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 &#119926;(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">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="keyword">const</span> <span class="special">&amp;</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">&lt;</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;&gt;</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">&lt;</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">)&gt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;</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">&lt;&lt;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;</span> <span class="keyword">const</span> <span class="special">&amp;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;::</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">&lt;</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">)&gt;</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">&lt;</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">)&gt;</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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;</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
      &#8719;<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 &#8719;<sub>i</sub> [ /a(1+&#949;)/, /b(1-&#949;)/ ]. (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 &#949; of the boundary is approximately &#949;<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">&lt;</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">)&gt;</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">&lt;</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">&gt;</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 &#169; 2006-2019 Nikhar
      Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
      Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
      R&#229;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="diff.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>