1 /*
2 * Copyright Nick Thompson, 2017
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 */
7
8 #include <boost/math/quadrature/naive_monte_carlo.hpp>
9 #include <iostream>
10 #include <iomanip>
11 #include <limits>
12 #include <cmath>
13 #include <thread>
14 #include <future>
15 #include <string>
16 #include <chrono>
17 #include <boost/math/special_functions/pow.hpp>
18 #include <boost/math/constants/constants.hpp>
19
20 using std::vector;
21 using std::pair;
22 using boost::math::quadrature::naive_monte_carlo;
23
display_progress(double progress,double error_estimate,double current_estimate,std::chrono::duration<double> estimated_time_to_completion)24 void display_progress(double progress,
25 double error_estimate,
26 double current_estimate,
27 std::chrono::duration<double> estimated_time_to_completion)
28 {
29 int barWidth = 70;
30
31 std::cout << "[";
32 int pos = barWidth * progress;
33 for (int i = 0; i < barWidth; ++i) {
34 if (i < pos) std::cout << "=";
35 else if (i == pos) std::cout << ">";
36 else std::cout << " ";
37 }
38 std::cout << "] "
39 << int(progress * 100.0)
40 << "%, E = "
41 << std::setprecision(3)
42 << error_estimate
43 << ", time to completion: "
44 << estimated_time_to_completion.count()
45 << " seconds, estimate: "
46 << std::setprecision(5)
47 << current_estimate
48 << " \r";
49
50 std::cout.flush();
51 }
52
main()53 int main()
54 {
55 double exact = 1.3932039296856768591842462603255;
56 double A = 1.0 / boost::math::pow<3>(boost::math::constants::pi<double>());
57 auto g = [&](std::vector<double> const & x)
58 {
59 return A / (1.0 - cos(x[0])*cos(x[1])*cos(x[2]));
60 };
61 vector<pair<double, double>> bounds{{0, boost::math::constants::pi<double>() }, {0, boost::math::constants::pi<double>() }, {0, boost::math::constants::pi<double>() }};
62 naive_monte_carlo<double, decltype(g)> mc(g, bounds, 0.001);
63
64 auto task = mc.integrate();
65
66 int s = 0;
67 std::cout << "Hit ctrl-c to cancel.\n";
68 while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready)
69 {
70 display_progress(mc.progress(),
71 mc.current_error_estimate(),
72 mc.current_estimate(),
73 mc.estimated_time_to_completion());
74 // TODO: The following shows that cancellation works,
75 // but it would be nice to show how it works with a ctrl-c signal handler.
76 if (s++ > 25){
77 mc.cancel();
78 std::cout << "\nCancelling because this is too slow!\n";
79 }
80 }
81 double y = task.get();
82 display_progress(mc.progress(),
83 mc.current_error_estimate(),
84 mc.current_estimate(),
85 mc.estimated_time_to_completion());
86 std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::fixed;
87 std::cout << "\nFinal value: " << y << std::endl;
88 std::cout << "Exact : " << exact << std::endl;
89 std::cout << "Final error estimate: " << mc.current_error_estimate() << std::endl;
90 std::cout << "Actual error : " << abs(y - exact) << std::endl;
91 std::cout << "Function calls: " << mc.calls() << std::endl;
92 std::cout << "Is this good enough? [y/N] ";
93 bool goodenough = true;
94 std::string line;
95 std::getline(std::cin, line);
96 if (line[0] != 'y')
97 {
98 goodenough = false;
99 }
100 double new_error = -1;
101 if (!goodenough)
102 {
103 std::cout << "What is the new target error? ";
104 std::getline(std::cin, line);
105 new_error = atof(line.c_str());
106 if (new_error >= mc.current_error_estimate())
107 {
108 std::cout << "That error bound is already satisfied.\n";
109 return 0;
110 }
111 }
112 if (new_error > 0)
113 {
114 mc.update_target_error(new_error);
115 auto task = mc.integrate();
116 std::cout << "Hit ctrl-c to cancel.\n";
117 while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready)
118 {
119 display_progress(mc.progress(),
120 mc.current_error_estimate(),
121 mc.current_estimate(),
122 mc.estimated_time_to_completion());
123 }
124 double y = task.get();
125 display_progress(mc.progress(),
126 mc.current_error_estimate(),
127 mc.current_estimate(),
128 mc.estimated_time_to_completion());
129 std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::fixed;
130 std::cout << "\nFinal value: " << y << std::endl;
131 std::cout << "Exact : " << exact << std::endl;
132 std::cout << "Final error estimate: " << mc.current_error_estimate() << std::endl;
133 std::cout << "Actual error : " << abs(y - exact) << std::endl;
134 std::cout << "Function calls: " << mc.calls() << std::endl;
135 }
136 }
137