• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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