469 lines · plain
1/*2 * Copyright Nick Thompson, 20183 * Use, modification and distribution are subject to the4 * Boost Software License, Version 1.0. (See accompanying file5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 */7#ifndef BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP8#define BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP9#include <sstream>10#include <algorithm>11#include <vector>12#include <atomic>13#include <memory>14#include <functional>15#include <future>16#include <thread>17#include <initializer_list>18#include <utility>19#include <random>20#include <chrono>21#include <map>22#include <type_traits>23#include <boost/math/policies/error_handling.hpp>24#include <boost/math/special_functions/fpclassify.hpp>25 26#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES27# include <iostream>28#endif29 30namespace boost { namespace math { namespace quadrature {31 32namespace detail {33 enum class limit_classification {FINITE,34 LOWER_BOUND_INFINITE,35 UPPER_BOUND_INFINITE,36 DOUBLE_INFINITE};37}38 39template<class Real, class F, class RandomNumberGenerator = std::mt19937_64, class Policy = boost::math::policies::policy<>,40 typename std::enable_if<std::is_trivially_copyable<Real>::value, bool>::type = true>41class naive_monte_carlo42{43public:44 naive_monte_carlo(const F& integrand,45 std::vector<std::pair<Real, Real>> const & bounds,46 Real error_goal,47 bool singular = true,48 uint64_t threads = std::thread::hardware_concurrency(),49 uint64_t seed = 0) noexcept : m_num_threads{threads}, m_seed{seed}, m_volume(1)50 {51 using std::numeric_limits;52 using std::sqrt;53 using boost::math::isinf;54 55 uint64_t n = bounds.size();56 m_lbs.resize(n);57 m_dxs.resize(n);58 m_limit_types.resize(n);59 60 static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";61 for (uint64_t i = 0; i < n; ++i)62 {63 if (bounds[i].second <= bounds[i].first)64 {65 boost::math::policies::raise_domain_error(function, "The upper bound is <= the lower bound.\n", bounds[i].second, Policy());66 return;67 }68 if (isinf(bounds[i].first))69 {70 if (isinf(bounds[i].second))71 {72 m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE;73 }74 else75 {76 m_limit_types[i] = detail::limit_classification::LOWER_BOUND_INFINITE;77 // Ok ok this is bad to use the second bound as the lower limit and then reflect.78 m_lbs[i] = bounds[i].second;79 m_dxs[i] = numeric_limits<Real>::quiet_NaN();80 }81 }82 else if (isinf(bounds[i].second))83 {84 m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE;85 if (singular)86 {87 // I've found that it's easier to sample on a closed set and perturb the boundary88 // than to try to sample very close to the boundary.89 m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)());90 }91 else92 {93 m_lbs[i] = bounds[i].first;94 }95 m_dxs[i] = numeric_limits<Real>::quiet_NaN();96 }97 else98 {99 m_limit_types[i] = detail::limit_classification::FINITE;100 if (singular)101 {102 if (bounds[i].first == 0)103 {104 m_lbs[i] = std::numeric_limits<Real>::epsilon();105 }106 else107 {108 m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)());109 }110 111 m_dxs[i] = std::nextafter(bounds[i].second, std::numeric_limits<Real>::lowest()) - m_lbs[i];112 }113 else114 {115 m_lbs[i] = bounds[i].first;116 m_dxs[i] = bounds[i].second - bounds[i].first;117 }118 m_volume *= m_dxs[i];119 }120 }121 122 m_integrand = [this, &integrand](std::vector<Real> & x)->Real123 {124 Real coeff = m_volume;125 for (uint64_t i = 0; i < x.size(); ++i)126 {127 // Variable transformation are listed at:128 // https://en.wikipedia.org/wiki/Numerical_integration129 // However, we've made some changes to these so that we can evaluate on a compact domain.130 if (m_limit_types[i] == detail::limit_classification::FINITE)131 {132 x[i] = m_lbs[i] + x[i]*m_dxs[i];133 }134 else if (m_limit_types[i] == detail::limit_classification::UPPER_BOUND_INFINITE)135 {136 Real t = x[i];137 Real z = 1/(1 + numeric_limits<Real>::epsilon() - t);138 coeff *= (z*z)*(1 + numeric_limits<Real>::epsilon());139 x[i] = m_lbs[i] + t*z;140 }141 else if (m_limit_types[i] == detail::limit_classification::LOWER_BOUND_INFINITE)142 {143 Real t = x[i];144 Real z = 1/(t+sqrt((numeric_limits<Real>::min)()));145 coeff *= (z*z);146 x[i] = m_lbs[i] + (t-1)*z;147 }148 else149 {150 Real t1 = 1/(1+numeric_limits<Real>::epsilon() - x[i]);151 Real t2 = 1/(x[i]+numeric_limits<Real>::epsilon());152 x[i] = (2*x[i]-1)*t1*t2/4;153 coeff *= (t1*t1+t2*t2)/4;154 }155 }156 return coeff*integrand(x);157 };158 159 // If we don't do a single function call in the constructor,160 // we can't do a restart.161 std::vector<Real> x(m_lbs.size());162 163 // If the seed is zero, that tells us to choose a random seed for the user:164 if (seed == 0)165 {166 std::random_device rd;167 seed = rd();168 }169 170 RandomNumberGenerator gen(seed);171 Real inv_denom = 1/static_cast<Real>(((gen.max)()-(gen.min)()));172 173 m_num_threads = (std::max)(m_num_threads, static_cast<uint64_t>(1));174 m_thread_calls.reset(new std::atomic<uint64_t>[threads]);175 m_thread_Ss.reset(new std::atomic<Real>[threads]);176 m_thread_averages.reset(new std::atomic<Real>[threads]);177 178 Real avg = 0;179 for (uint64_t i = 0; i < m_num_threads; ++i)180 {181 for (uint64_t j = 0; j < m_lbs.size(); ++j)182 {183 x[j] = (gen()-(gen.min)())*inv_denom;184 }185 Real y = m_integrand(x);186 m_thread_averages[i] = y; // relaxed store187 m_thread_calls[i] = 1;188 m_thread_Ss[i] = 0;189 avg += y;190 }191 avg /= m_num_threads;192 m_avg = avg; // relaxed store193 194 m_error_goal = error_goal; // relaxed store195 m_start = std::chrono::system_clock::now();196 m_done = false; // relaxed store197 m_total_calls = m_num_threads; // relaxed store198 m_variance = (numeric_limits<Real>::max)();199 }200 201 std::future<Real> integrate()202 {203 // Set done to false in case we wish to restart:204 m_done.store(false); // relaxed store, no worker threads yet205 m_start = std::chrono::system_clock::now();206 return std::async(std::launch::async,207 &naive_monte_carlo::m_integrate, this);208 }209 210 void cancel()211 {212 // If seed = 0 (meaning have the routine pick the seed), this leaves the seed the same.213 // If seed != 0, then the seed is changed, so a restart doesn't do the exact same thing.214 m_seed = m_seed*m_seed;215 m_done = true; // relaxed store, worker threads will get the message eventually216 // Make sure the error goal is infinite, because otherwise we'll loop when we do the final error goal check:217 m_error_goal = (std::numeric_limits<Real>::max)();218 }219 220 Real variance() const221 {222 return m_variance.load();223 }224 225 Real current_error_estimate() const226 {227 using std::sqrt;228 //229 // There is a bug here: m_variance and m_total_calls get updated asynchronously230 // and may be out of synch when we compute the error estimate, not sure if it matters though...231 //232 return sqrt(m_variance.load()/m_total_calls.load());233 }234 235 std::chrono::duration<Real> estimated_time_to_completion() const236 {237 auto now = std::chrono::system_clock::now();238 std::chrono::duration<Real> elapsed_seconds = now - m_start;239 Real r = this->current_error_estimate()/m_error_goal.load(); // relaxed load240 if (r*r <= 1) {241 return 0*elapsed_seconds;242 }243 return (r*r - 1)*elapsed_seconds;244 }245 246 void update_target_error(Real new_target_error)247 {248 m_error_goal = new_target_error; // relaxed store249 }250 251 Real progress() const252 {253 Real r = m_error_goal.load()/this->current_error_estimate(); // relaxed load254 if (r*r >= 1)255 {256 return 1;257 }258 return r*r;259 }260 261 Real current_estimate() const262 {263 return m_avg.load();264 }265 266 uint64_t calls() const267 {268 return m_total_calls.load(); // relaxed load269 }270 271private:272 273 Real m_integrate()274 {275 uint64_t seed;276 // If the user tells us to pick a seed, pick a seed:277 if (m_seed == 0)278 {279 std::random_device rd;280 seed = rd();281 }282 else // use the seed we are given:283 {284 seed = m_seed;285 }286 RandomNumberGenerator gen(seed);287 int max_repeat_tries = 5;288 do{289 290 if (max_repeat_tries < 5)291 {292 m_done = false;293 294#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES295 std::cerr << "Failed to achieve required tolerance first time through..\n";296 std::cerr << " variance = " << m_variance << std::endl;297 std::cerr << " average = " << m_avg << std::endl;298 std::cerr << " total calls = " << m_total_calls << std::endl;299 300 for (std::size_t i = 0; i < m_num_threads; ++i)301 std::cerr << " thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;302 for (std::size_t i = 0; i < m_num_threads; ++i)303 std::cerr << " thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;304 for (std::size_t i = 0; i < m_num_threads; ++i)305 std::cerr << " thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;306#endif307 }308 309 std::vector<std::thread> threads(m_num_threads);310 for (uint64_t i = 0; i < threads.size(); ++i)311 {312 threads[i] = std::thread(&naive_monte_carlo::m_thread_monte, this, i, gen());313 }314 do {315 std::this_thread::sleep_for(std::chrono::milliseconds(100));316 uint64_t total_calls = 0;317 for (uint64_t i = 0; i < m_num_threads; ++i)318 {319 uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);320 total_calls += t_calls;321 }322 Real variance = 0;323 Real avg = 0;324 for (uint64_t i = 0; i < m_num_threads; ++i)325 {326 uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);327 // Will this overflow? Not hard to remove . . .328 avg += m_thread_averages[i].load(std::memory_order_relaxed)*(static_cast<Real>(t_calls) / static_cast<Real>(total_calls));329 variance += m_thread_Ss[i].load(std::memory_order_relaxed);330 }331 m_avg.store(avg, std::memory_order_release);332 m_variance.store(variance / (total_calls - 1), std::memory_order_release);333 m_total_calls = total_calls; // relaxed store, it's just for user feedback334 // Allow cancellation:335 if (m_done) // relaxed load336 {337 break;338 }339 } while (m_total_calls < 2048 || this->current_error_estimate() > m_error_goal.load(std::memory_order_consume));340 // Error bound met; signal the threads:341 m_done = true; // relaxed store, threads will get the message in the end342 std::for_each(threads.begin(), threads.end(),343 std::mem_fn(&std::thread::join));344 if (m_exception)345 {346 std::rethrow_exception(m_exception);347 }348 // Incorporate their work into the final estimate:349 uint64_t total_calls = 0;350 for (uint64_t i = 0; i < m_num_threads; ++i)351 {352 uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);353 total_calls += t_calls;354 }355 Real variance = 0;356 Real avg = 0;357 358 for (uint64_t i = 0; i < m_num_threads; ++i)359 {360 uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);361 // Averages weighted by the number of calls the thread made:362 avg += m_thread_averages[i].load(std::memory_order_relaxed)*(static_cast<Real>(t_calls) / static_cast<Real>(total_calls));363 variance += m_thread_Ss[i].load(std::memory_order_relaxed);364 }365 m_avg.store(avg, std::memory_order_release);366 m_variance.store(variance / (total_calls - 1), std::memory_order_release);367 m_total_calls = total_calls; // relaxed store, this is just user feedback368 369 // Sometimes, the master will observe the variance at a very "good" (or bad?) moment,370 // Then the threads proceed to find the variance is much greater by the time they hear the message to stop.371 // This *WOULD* make sure that the final error estimate is within the error bounds.372 }373 while ((--max_repeat_tries >= 0) && (this->current_error_estimate() > m_error_goal));374 375 return m_avg.load(std::memory_order_consume);376 }377 378 void m_thread_monte(uint64_t thread_index, uint64_t seed)379 {380 using std::numeric_limits;381 try382 {383 std::vector<Real> x(m_lbs.size());384 RandomNumberGenerator gen(seed);385 Real inv_denom = static_cast<Real>(1) / static_cast<Real>(( (gen.max)() - (gen.min)() ));386 Real M1 = m_thread_averages[thread_index].load(std::memory_order_consume);387 Real S = m_thread_Ss[thread_index].load(std::memory_order_consume);388 // Kahan summation is required or the value of the integrand will go on a random walk during long computations.389 // See the implementation discussion.390 // The idea is that the unstabilized additions have error sigma(f)/sqrt(N) + epsilon*N, which diverges faster than it converges!391 // Kahan summation turns this to sigma(f)/sqrt(N) + epsilon^2*N, and the random walk occurs on a timescale of 10^14 years (on current hardware)392 Real compensator = 0;393 uint64_t k = m_thread_calls[thread_index].load(std::memory_order_consume);394 while (!m_done) // relaxed load395 {396 int j = 0;397 // If we don't have a certain number of calls before an update, we can easily terminate prematurely398 // because the variance estimate is way too low. This magic number is a reasonable compromise, as 1/sqrt(2048) = 0.02,399 // so it should recover 2 digits if the integrand isn't poorly behaved, and if it is, it should discover that before premature termination.400 // Of course if the user has 64 threads, then this number is probably excessive.401 int magic_calls_before_update = 2048;402 while (j++ < magic_calls_before_update)403 {404 for (uint64_t i = 0; i < m_lbs.size(); ++i)405 {406 x[i] = (gen() - (gen.min)())*inv_denom;407 }408 Real f = m_integrand(x);409 using std::isfinite;410 if (!isfinite(f))411 {412 // The call to m_integrand transform x, so this error message states the correct node.413 std::stringstream os;414 os << "Your integrand was evaluated at {";415 for (uint64_t i = 0; i < x.size() -1; ++i)416 {417 os << x[i] << ", ";418 }419 os << x[x.size() -1] << "}, and returned " << f << std::endl;420 static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";421 boost::math::policies::raise_domain_error(function, os.str().c_str(), /*this is a dummy arg to make it compile*/ 7.2, Policy());422 }423 ++k;424 Real term = (f - M1)/k;425 Real y1 = term - compensator;426 Real M2 = M1 + y1;427 compensator = (M2 - M1) - y1;428 S += (f - M1)*(f - M2);429 M1 = M2;430 }431 m_thread_averages[thread_index].store(M1, std::memory_order_release);432 m_thread_Ss[thread_index].store(S, std::memory_order_release);433 m_thread_calls[thread_index].store(k, std::memory_order_release);434 }435 }436 catch (...)437 {438 // Signal the other threads that the computation is ruined:439 m_done = true; // relaxed store440 std::lock_guard<std::mutex> lock(m_exception_mutex); // Scoped lock to prevent race writing to m_exception441 m_exception = std::current_exception();442 }443 }444 445 std::function<Real(std::vector<Real> &)> m_integrand;446 uint64_t m_num_threads;447 std::atomic<uint64_t> m_seed;448 std::atomic<Real> m_error_goal;449 std::atomic<bool> m_done{};450 std::vector<Real> m_lbs;451 std::vector<Real> m_dxs;452 std::vector<detail::limit_classification> m_limit_types;453 Real m_volume;454 std::atomic<uint64_t> m_total_calls{};455 // I wanted these to be vectors rather than maps,456 // but you can't resize a vector of atomics.457 std::unique_ptr<std::atomic<uint64_t>[]> m_thread_calls;458 std::atomic<Real> m_variance;459 std::unique_ptr<std::atomic<Real>[]> m_thread_Ss;460 std::atomic<Real> m_avg;461 std::unique_ptr<std::atomic<Real>[]> m_thread_averages;462 std::chrono::time_point<std::chrono::system_clock> m_start;463 std::exception_ptr m_exception;464 std::mutex m_exception_mutex;465};466 467}}}468#endif469