123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136 |
- /*
- * Copyright Nick Thompson, 2017
- * Use, modification and distribution are subject to the
- * Boost Software License, Version 1.0. (See accompanying file
- * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- */
- #include <boost/math/quadrature/naive_monte_carlo.hpp>
- #include <iostream>
- #include <iomanip>
- #include <limits>
- #include <cmath>
- #include <thread>
- #include <future>
- #include <string>
- #include <chrono>
- #include <boost/math/special_functions/pow.hpp>
- #include <boost/math/constants/constants.hpp>
- using std::vector;
- using std::pair;
- using boost::math::quadrature::naive_monte_carlo;
- void display_progress(double progress,
- double error_estimate,
- double current_estimate,
- std::chrono::duration<double> estimated_time_to_completion)
- {
- int barWidth = 70;
- std::cout << "[";
- int pos = barWidth * progress;
- for (int i = 0; i < barWidth; ++i) {
- if (i < pos) std::cout << "=";
- else if (i == pos) std::cout << ">";
- else std::cout << " ";
- }
- std::cout << "] "
- << int(progress * 100.0)
- << "%, E = "
- << std::setprecision(3)
- << error_estimate
- << ", time to completion: "
- << estimated_time_to_completion.count()
- << " seconds, estimate: "
- << std::setprecision(5)
- << current_estimate
- << " \r";
- std::cout.flush();
- }
- int main()
- {
- double exact = 1.3932039296856768591842462603255;
- double A = 1.0 / boost::math::pow<3>(boost::math::constants::pi<double>());
- auto g = [&](std::vector<double> const & x)
- {
- return A / (1.0 - cos(x[0])*cos(x[1])*cos(x[2]));
- };
- vector<pair<double, double>> bounds{{0, boost::math::constants::pi<double>() }, {0, boost::math::constants::pi<double>() }, {0, boost::math::constants::pi<double>() }};
- naive_monte_carlo<double, decltype(g)> mc(g, bounds, 0.001);
- auto task = mc.integrate();
- int s = 0;
- std::cout << "Hit ctrl-c to cancel.\n";
- while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready)
- {
- display_progress(mc.progress(),
- mc.current_error_estimate(),
- mc.current_estimate(),
- mc.estimated_time_to_completion());
- // TODO: The following shows that cancellation works,
- // but it would be nice to show how it works with a ctrl-c signal handler.
- if (s++ > 25){
- mc.cancel();
- std::cout << "\nCancelling because this is too slow!\n";
- }
- }
- double y = task.get();
- display_progress(mc.progress(),
- mc.current_error_estimate(),
- mc.current_estimate(),
- mc.estimated_time_to_completion());
- std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::fixed;
- std::cout << "\nFinal value: " << y << std::endl;
- std::cout << "Exact : " << exact << std::endl;
- std::cout << "Final error estimate: " << mc.current_error_estimate() << std::endl;
- std::cout << "Actual error : " << abs(y - exact) << std::endl;
- std::cout << "Function calls: " << mc.calls() << std::endl;
- std::cout << "Is this good enough? [y/N] ";
- bool goodenough = true;
- std::string line;
- std::getline(std::cin, line);
- if (line[0] != 'y')
- {
- goodenough = false;
- }
- double new_error = -1;
- if (!goodenough)
- {
- std::cout << "What is the new target error? ";
- std::getline(std::cin, line);
- new_error = atof(line.c_str());
- if (new_error >= mc.current_error_estimate())
- {
- std::cout << "That error bound is already satisfied.\n";
- return 0;
- }
- }
- if (new_error > 0)
- {
- mc.update_target_error(new_error);
- auto task = mc.integrate();
- std::cout << "Hit ctrl-c to cancel.\n";
- while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready)
- {
- display_progress(mc.progress(),
- mc.current_error_estimate(),
- mc.current_estimate(),
- mc.estimated_time_to_completion());
- }
- double y = task.get();
- display_progress(mc.progress(),
- mc.current_error_estimate(),
- mc.current_estimate(),
- mc.estimated_time_to_completion());
- std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::fixed;
- std::cout << "\nFinal value: " << y << std::endl;
- std::cout << "Exact : " << exact << std::endl;
- std::cout << "Final error estimate: " << mc.current_error_estimate() << std::endl;
- std::cout << "Actual error : " << abs(y - exact) << std::endl;
- std::cout << "Function calls: " << mc.calls() << std::endl;
- }
- }
|