]>
Commit | Line | Data |
---|---|---|
20effc67 TL |
1 | // (C) Copyright Nick Thompson 2020. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | #ifndef BOOST_MATH_TOOLS_ULP_PLOT_HPP | |
6 | #define BOOST_MATH_TOOLS_ULP_PLOT_HPP | |
7 | #include <algorithm> | |
8 | #include <iostream> | |
9 | #include <iomanip> | |
10 | #include <cassert> | |
11 | #include <vector> | |
12 | #include <utility> | |
13 | #include <fstream> | |
14 | #include <string> | |
15 | #include <list> | |
16 | #include <random> | |
17 | #include <stdexcept> | |
18 | #include <boost/math/tools/condition_numbers.hpp> | |
19 | #include <boost/random/uniform_real_distribution.hpp> | |
20 | #include <boost/algorithm/string/predicate.hpp> | |
21 | ||
22 | ||
23 | // Design of this function comes from: | |
24 | // https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/ | |
25 | ||
26 | // The envelope is the maximum of 1/2 and half the condition number of function evaluation. | |
27 | ||
28 | namespace boost::math::tools { | |
29 | ||
30 | namespace detail { | |
31 | template<class F1, class F2, class CoarseReal, class PreciseReal> | |
32 | void write_gridlines(std::ostream& fs, int horizontal_lines, int vertical_lines, | |
33 | F1 x_scale, F2 y_scale, CoarseReal min_x, CoarseReal max_x, PreciseReal min_y, PreciseReal max_y, | |
34 | int graph_width, int graph_height, int margin_left, std::string const & font_color) | |
35 | { | |
36 | // Make a grid: | |
37 | for (int i = 1; i <= horizontal_lines; ++i) { | |
38 | PreciseReal y_cord_dataspace = min_y + ((max_y - min_y)*i)/horizontal_lines; | |
39 | auto y = y_scale(y_cord_dataspace); | |
40 | fs << "<line x1='0' y1='" << y << "' x2='" << graph_width | |
41 | << "' y2='" << y | |
42 | << "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n"; | |
43 | ||
44 | fs << "<text x='" << -margin_left/4 + 5 << "' y='" << y - 3 | |
45 | << "' font-family='times' font-size='10' fill='" << font_color << "' transform='rotate(-90 " | |
46 | << -margin_left/4 + 8 << " " << y + 5 << ")'>" | |
47 | << std::setprecision(4) << y_cord_dataspace << "</text>\n"; | |
48 | } | |
49 | ||
50 | for (int i = 1; i <= vertical_lines; ++i) { | |
51 | CoarseReal x_cord_dataspace = min_x + ((max_x - min_x)*i)/vertical_lines; | |
52 | CoarseReal x = x_scale(x_cord_dataspace); | |
53 | fs << "<line x1='" << x << "' y1='0' x2='" << x | |
54 | << "' y2='" << graph_height | |
55 | << "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n"; | |
56 | ||
57 | fs << "<text x='" << x - 10 << "' y='" << graph_height + 10 | |
58 | << "' font-family='times' font-size='10' fill='" << font_color << "'>" | |
59 | << std::setprecision(4) << x_cord_dataspace << "</text>\n"; | |
60 | } | |
61 | } | |
62 | } | |
63 | ||
64 | template<class F, typename PreciseReal, typename CoarseReal> | |
65 | class ulps_plot { | |
66 | public: | |
67 | ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, | |
68 | size_t samples = 1000, bool perturb_abscissas = false, int random_seed = -1); | |
69 | ||
70 | ulps_plot& clip(PreciseReal clip); | |
71 | ||
72 | ulps_plot& width(int width); | |
73 | ||
74 | ulps_plot& envelope_color(std::string const & color); | |
75 | ||
76 | ulps_plot& title(std::string const & title); | |
77 | ||
78 | ulps_plot& background_color(std::string const & background_color); | |
79 | ||
80 | ulps_plot& font_color(std::string const & font_color); | |
81 | ||
82 | ulps_plot& crop_color(std::string const & color); | |
83 | ||
84 | ulps_plot& nan_color(std::string const & color); | |
85 | ||
86 | ulps_plot& ulp_envelope(bool write_ulp); | |
87 | ||
88 | template<class G> | |
89 | ulps_plot& add_fn(G g, std::string const & color = "steelblue"); | |
90 | ||
91 | ulps_plot& horizontal_lines(int horizontal_lines); | |
92 | ||
93 | ulps_plot& vertical_lines(int vertical_lines); | |
94 | ||
95 | void write(std::string const & filename) const; | |
96 | ||
97 | friend std::ostream& operator<<(std::ostream& fs, ulps_plot const & plot) | |
98 | { | |
99 | using std::abs; | |
100 | using std::floor; | |
101 | using std::isnan; | |
102 | if (plot.ulp_list_.size() == 0) | |
103 | { | |
104 | throw std::domain_error("No functions added for comparison."); | |
105 | } | |
106 | if (plot.width_ <= 1) | |
107 | { | |
108 | throw std::domain_error("Width = " + std::to_string(plot.width_) + ", which is too small."); | |
109 | } | |
110 | ||
111 | PreciseReal worst_ulp_distance = 0; | |
112 | PreciseReal min_y = std::numeric_limits<PreciseReal>::max(); | |
113 | PreciseReal max_y = std::numeric_limits<PreciseReal>::lowest(); | |
114 | for (auto const & ulp_vec : plot.ulp_list_) | |
115 | { | |
116 | for (auto const & ulp : ulp_vec) | |
117 | { | |
118 | if (static_cast<PreciseReal>(abs(ulp)) > worst_ulp_distance) | |
119 | { | |
120 | worst_ulp_distance = static_cast<PreciseReal>(abs(ulp)); | |
121 | } | |
122 | if (static_cast<PreciseReal>(ulp) < min_y) | |
123 | { | |
124 | min_y = static_cast<PreciseReal>(ulp); | |
125 | } | |
126 | if (static_cast<PreciseReal>(ulp) > max_y) | |
127 | { | |
128 | max_y = static_cast<PreciseReal>(ulp); | |
129 | } | |
130 | } | |
131 | } | |
132 | ||
133 | // half-ulp accuracy is the best that can be expected; sometimes we can get less, but barely less. | |
134 | // then the axes don't show up; painful! | |
135 | if (max_y < 0.5) { | |
136 | max_y = 0.5; | |
137 | } | |
138 | if (min_y > -0.5) { | |
139 | min_y = -0.5; | |
140 | } | |
141 | ||
142 | if (plot.clip_ > 0) | |
143 | { | |
144 | if (max_y > plot.clip_) | |
145 | { | |
146 | max_y = plot.clip_; | |
147 | } | |
148 | if (min_y < -plot.clip_) | |
149 | { | |
150 | min_y = -plot.clip_; | |
151 | } | |
152 | } | |
153 | ||
154 | int height = static_cast<int>(floor(double(plot.width_)/1.61803)); | |
155 | int margin_top = 40; | |
156 | int margin_left = 25; | |
157 | if (plot.title_.size() == 0) | |
158 | { | |
159 | margin_top = 10; | |
160 | margin_left = 15; | |
161 | } | |
162 | int margin_bottom = 20; | |
163 | int margin_right = 20; | |
164 | int graph_height = height - margin_bottom - margin_top; | |
165 | int graph_width = plot.width_ - margin_left - margin_right; | |
166 | ||
167 | // Maps [a,b] to [0, graph_width] | |
168 | auto x_scale = [&](CoarseReal x)->CoarseReal | |
169 | { | |
170 | return ((x-plot.a_)/(plot.b_ - plot.a_))*static_cast<CoarseReal>(graph_width); | |
171 | }; | |
172 | ||
173 | auto y_scale = [&](PreciseReal y)->PreciseReal | |
174 | { | |
175 | return ((max_y - y)/(max_y - min_y) )*static_cast<PreciseReal>(graph_height); | |
176 | }; | |
177 | ||
178 | fs << "<?xml version=\"1.0\" encoding='UTF-8' ?>\n" | |
179 | << "<svg xmlns='http://www.w3.org/2000/svg' width='" | |
180 | << plot.width_ << "' height='" | |
181 | << height << "'>\n" | |
182 | << "<style>\nsvg { background-color:" << plot.background_color_ << "; }\n" | |
183 | << "</style>\n"; | |
184 | if (plot.title_.size() > 0) | |
185 | { | |
186 | fs << "<text x='" << floor(plot.width_/2) | |
187 | << "' y='" << floor(margin_top/2) | |
188 | << "' font-family='Palatino' font-size='25' fill='" | |
189 | << plot.font_color_ << "' alignment-baseline='middle' text-anchor='middle'>" | |
190 | << plot.title_ | |
191 | << "</text>\n"; | |
192 | } | |
193 | ||
194 | // Construct SVG group to simplify the calculations slightly: | |
195 | fs << "<g transform='translate(" << margin_left << ", " << margin_top << ")'>\n"; | |
196 | // y-axis: | |
197 | fs << "<line x1='0' y1='0' x2='0' y2='" << graph_height | |
198 | << "' stroke='gray' stroke-width='1'/>\n"; | |
199 | PreciseReal x_axis_loc = y_scale(static_cast<PreciseReal>(0)); | |
200 | fs << "<line x1='0' y1='" << x_axis_loc | |
201 | << "' x2='" << graph_width << "' y2='" << x_axis_loc | |
202 | << "' stroke='gray' stroke-width='1'/>\n"; | |
203 | ||
204 | if (worst_ulp_distance > 3) | |
205 | { | |
206 | detail::write_gridlines(fs, plot.horizontal_lines_, plot.vertical_lines_, x_scale, y_scale, plot.a_, plot.b_, | |
207 | min_y, max_y, graph_width, graph_height, margin_left, plot.font_color_); | |
208 | } | |
209 | else | |
210 | { | |
211 | std::vector<double> ys{-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}; | |
212 | for (size_t i = 0; i < ys.size(); ++i) | |
213 | { | |
214 | if (min_y <= ys[i] && ys[i] <= max_y) | |
215 | { | |
216 | PreciseReal y_cord_dataspace = ys[i]; | |
217 | PreciseReal y = y_scale(y_cord_dataspace); | |
218 | fs << "<line x1='0' y1='" << y << "' x2='" << graph_width | |
219 | << "' y2='" << y | |
220 | << "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n"; | |
221 | ||
222 | fs << "<text x='" << -margin_left/2 << "' y='" << y - 3 | |
223 | << "' font-family='times' font-size='10' fill='" << plot.font_color_ << "' transform='rotate(-90 " | |
224 | << -margin_left/2 + 7 << " " << y << ")'>" | |
225 | << std::setprecision(4) << y_cord_dataspace << "</text>\n"; | |
226 | } | |
227 | } | |
228 | for (int i = 1; i <= plot.vertical_lines_; ++i) | |
229 | { | |
230 | CoarseReal x_cord_dataspace = plot.a_ + ((plot.b_ - plot.a_)*i)/plot.vertical_lines_; | |
231 | CoarseReal x = x_scale(x_cord_dataspace); | |
232 | fs << "<line x1='" << x << "' y1='0' x2='" << x | |
233 | << "' y2='" << graph_height | |
234 | << "' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />\n"; | |
235 | ||
236 | fs << "<text x='" << x - 10 << "' y='" << graph_height + 10 | |
237 | << "' font-family='times' font-size='10' fill='" << plot.font_color_ << "'>" | |
238 | << std::setprecision(4) << x_cord_dataspace << "</text>\n"; | |
239 | } | |
240 | } | |
241 | ||
242 | int color_idx = 0; | |
243 | for (auto const & ulp : plot.ulp_list_) | |
244 | { | |
245 | std::string color = plot.colors_[color_idx++]; | |
246 | for (size_t j = 0; j < ulp.size(); ++j) | |
247 | { | |
248 | if (isnan(ulp[j])) | |
249 | { | |
250 | if(plot.nan_color_ == "") | |
251 | continue; | |
252 | CoarseReal x = x_scale(plot.coarse_abscissas_[j]); | |
253 | PreciseReal y = y_scale(static_cast<PreciseReal>(plot.clip_)); | |
254 | fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.nan_color_ << "'/>\n"; | |
255 | y = y_scale(static_cast<PreciseReal>(-plot.clip_)); | |
256 | fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.nan_color_ << "'/>\n"; | |
257 | } | |
258 | if (plot.clip_ > 0 && static_cast<PreciseReal>(abs(ulp[j])) > plot.clip_) | |
259 | { | |
260 | if (plot.crop_color_ == "") | |
261 | continue; | |
262 | CoarseReal x = x_scale(plot.coarse_abscissas_[j]); | |
263 | PreciseReal y = y_scale(static_cast<PreciseReal>(ulp[j] < 0 ? -plot.clip_ : plot.clip_)); | |
264 | fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.crop_color_ << "'/>\n"; | |
265 | } | |
266 | else | |
267 | { | |
268 | CoarseReal x = x_scale(plot.coarse_abscissas_[j]); | |
269 | PreciseReal y = y_scale(static_cast<PreciseReal>(ulp[j])); | |
270 | fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << color << "'/>\n"; | |
271 | } | |
272 | } | |
273 | } | |
274 | ||
275 | if (plot.ulp_envelope_) | |
276 | { | |
277 | std::string close_path = "' stroke='" + plot.envelope_color_ + "' stroke-width='1' fill='none'></path>\n"; | |
278 | size_t jstart = 0; | |
279 | while (plot.cond_[jstart] > max_y) | |
280 | { | |
281 | ++jstart; | |
282 | if (jstart >= plot.cond_.size()) | |
283 | { | |
284 | goto done; | |
285 | } | |
286 | } | |
287 | ||
288 | size_t jmin = jstart; | |
289 | new_top_path: | |
290 | if (jmin >= plot.cond_.size()) | |
291 | { | |
292 | goto start_bottom_paths; | |
293 | } | |
294 | fs << "<path d='M" << x_scale(plot.coarse_abscissas_[jmin]) << " " << y_scale(plot.cond_[jmin]); | |
295 | ||
296 | for (size_t j = jmin + 1; j < plot.coarse_abscissas_.size(); ++j) | |
297 | { | |
298 | bool bad = isnan(plot.cond_[j]) || (plot.cond_[j] > max_y); | |
299 | if (bad) | |
300 | { | |
301 | ++j; | |
302 | while ( (j < plot.coarse_abscissas_.size() - 2) && bad) | |
303 | { | |
304 | bad = isnan(plot.cond_[j]) || (plot.cond_[j] > max_y); | |
305 | ++j; | |
306 | } | |
307 | jmin = j; | |
308 | fs << close_path; | |
309 | goto new_top_path; | |
310 | } | |
311 | ||
312 | CoarseReal t = x_scale(plot.coarse_abscissas_[j]); | |
313 | PreciseReal y = y_scale(plot.cond_[j]); | |
314 | fs << " L" << t << " " << y; | |
315 | } | |
316 | fs << close_path; | |
317 | start_bottom_paths: | |
318 | jmin = jstart; | |
319 | new_bottom_path: | |
320 | if (jmin >= plot.cond_.size()) | |
321 | { | |
322 | goto done; | |
323 | } | |
324 | fs << "<path d='M" << x_scale(plot.coarse_abscissas_[jmin]) << " " << y_scale(-plot.cond_[jmin]); | |
325 | ||
326 | for (size_t j = jmin + 1; j < plot.coarse_abscissas_.size(); ++j) | |
327 | { | |
328 | bool bad = isnan(plot.cond_[j]) || (-plot.cond_[j] < min_y); | |
329 | if (bad) | |
330 | { | |
331 | ++j; | |
332 | while ( (j < plot.coarse_abscissas_.size() - 2) && bad) | |
333 | { | |
334 | bad = isnan(plot.cond_[j]) || (-plot.cond_[j] < min_y); | |
335 | ++j; | |
336 | } | |
337 | jmin = j; | |
338 | fs << close_path; | |
339 | goto new_bottom_path; | |
340 | } | |
341 | CoarseReal t = x_scale(plot.coarse_abscissas_[j]); | |
342 | PreciseReal y = y_scale(-plot.cond_[j]); | |
343 | fs << " L" << t << " " << y; | |
344 | } | |
345 | fs << close_path; | |
346 | } | |
347 | done: | |
348 | fs << "</g>\n" | |
349 | << "</svg>\n"; | |
350 | return fs; | |
351 | } | |
352 | ||
353 | private: | |
354 | std::vector<PreciseReal> precise_abscissas_; | |
355 | std::vector<CoarseReal> coarse_abscissas_; | |
356 | std::vector<PreciseReal> precise_ordinates_; | |
357 | std::vector<PreciseReal> cond_; | |
358 | std::list<std::vector<CoarseReal>> ulp_list_; | |
359 | std::vector<std::string> colors_; | |
360 | CoarseReal a_; | |
361 | CoarseReal b_; | |
362 | PreciseReal clip_; | |
363 | int width_; | |
364 | std::string envelope_color_; | |
365 | bool ulp_envelope_; | |
366 | int horizontal_lines_; | |
367 | int vertical_lines_; | |
368 | std::string title_; | |
369 | std::string background_color_; | |
370 | std::string font_color_; | |
371 | std::string crop_color_; | |
372 | std::string nan_color_; | |
373 | }; | |
374 | ||
375 | template<class F, typename PreciseReal, typename CoarseReal> | |
376 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::envelope_color(std::string const & color) | |
377 | { | |
378 | envelope_color_ = color; | |
379 | return *this; | |
380 | } | |
381 | ||
382 | template<class F, typename PreciseReal, typename CoarseReal> | |
383 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::clip(PreciseReal clip) | |
384 | { | |
385 | clip_ = clip; | |
386 | return *this; | |
387 | } | |
388 | ||
389 | template<class F, typename PreciseReal, typename CoarseReal> | |
390 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::width(int width) | |
391 | { | |
392 | width_ = width; | |
393 | return *this; | |
394 | } | |
395 | ||
396 | template<class F, typename PreciseReal, typename CoarseReal> | |
397 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::horizontal_lines(int horizontal_lines) | |
398 | { | |
399 | horizontal_lines_ = horizontal_lines; | |
400 | return *this; | |
401 | } | |
402 | ||
403 | template<class F, typename PreciseReal, typename CoarseReal> | |
404 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::vertical_lines(int vertical_lines) | |
405 | { | |
406 | vertical_lines_ = vertical_lines; | |
407 | return *this; | |
408 | } | |
409 | ||
410 | template<class F, typename PreciseReal, typename CoarseReal> | |
411 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::title(std::string const & title) | |
412 | { | |
413 | title_ = title; | |
414 | return *this; | |
415 | } | |
416 | ||
417 | template<class F, typename PreciseReal, typename CoarseReal> | |
418 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::background_color(std::string const & background_color) | |
419 | { | |
420 | background_color_ = background_color; | |
421 | return *this; | |
422 | } | |
423 | ||
424 | template<class F, typename PreciseReal, typename CoarseReal> | |
425 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::font_color(std::string const & font_color) | |
426 | { | |
427 | font_color_ = font_color; | |
428 | return *this; | |
429 | } | |
430 | ||
431 | template<class F, typename PreciseReal, typename CoarseReal> | |
432 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::crop_color(std::string const & color) | |
433 | { | |
434 | crop_color_ = color; | |
435 | return *this; | |
436 | } | |
437 | ||
438 | template<class F, typename PreciseReal, typename CoarseReal> | |
439 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::nan_color(std::string const & color) | |
440 | { | |
441 | nan_color_ = color; | |
442 | return *this; | |
443 | } | |
444 | ||
445 | template<class F, typename PreciseReal, typename CoarseReal> | |
446 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::ulp_envelope(bool write_ulp_envelope) | |
447 | { | |
448 | ulp_envelope_ = write_ulp_envelope; | |
449 | return *this; | |
450 | } | |
451 | ||
452 | template<class F, typename PreciseReal, typename CoarseReal> | |
453 | void ulps_plot<F, PreciseReal, CoarseReal>::write(std::string const & filename) const | |
454 | { | |
455 | if (!boost::algorithm::ends_with(filename, ".svg")) | |
456 | { | |
457 | throw std::logic_error("Only svg files are supported at this time."); | |
458 | } | |
459 | std::ofstream fs(filename); | |
460 | fs << *this; | |
461 | fs.close(); | |
462 | } | |
463 | ||
464 | ||
465 | template<class F, typename PreciseReal, typename CoarseReal> | |
466 | ulps_plot<F, PreciseReal, CoarseReal>::ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, | |
467 | size_t samples, bool perturb_abscissas, int random_seed) : crop_color_("red") | |
468 | { | |
469 | // Use digits10 for this comparison in case the two types have differeing radixes: | |
470 | static_assert(std::numeric_limits<PreciseReal>::digits10 >= std::numeric_limits<CoarseReal>::digits10, "PreciseReal must have higher precision that CoarseReal"); | |
471 | if (samples < 10) | |
472 | { | |
473 | throw std::domain_error("Must have at least 10 samples, samples = " + std::to_string(samples)); | |
474 | } | |
475 | if (b <= a) | |
476 | { | |
477 | throw std::domain_error("On interval [a,b], b > a is required."); | |
478 | } | |
479 | a_ = a; | |
480 | b_ = b; | |
481 | ||
482 | std::mt19937_64 gen; | |
483 | if (random_seed == -1) | |
484 | { | |
485 | std::random_device rd; | |
486 | gen.seed(rd()); | |
487 | } | |
488 | // Boost's uniform_real_distribution can generate quad and multiprecision random numbers; std's cannot: | |
489 | boost::random::uniform_real_distribution<PreciseReal> dis(static_cast<PreciseReal>(a), static_cast<PreciseReal>(b)); | |
490 | precise_abscissas_.resize(samples); | |
491 | coarse_abscissas_.resize(samples); | |
492 | ||
493 | if (perturb_abscissas) | |
494 | { | |
495 | for(size_t i = 0; i < samples; ++i) | |
496 | { | |
497 | precise_abscissas_[i] = dis(gen); | |
498 | } | |
499 | std::sort(precise_abscissas_.begin(), precise_abscissas_.end()); | |
500 | for (size_t i = 0; i < samples; ++i) | |
501 | { | |
502 | coarse_abscissas_[i] = static_cast<CoarseReal>(precise_abscissas_[i]); | |
503 | } | |
504 | } | |
505 | else | |
506 | { | |
507 | for(size_t i = 0; i < samples; ++i) | |
508 | { | |
509 | coarse_abscissas_[i] = static_cast<CoarseReal>(dis(gen)); | |
510 | } | |
511 | std::sort(coarse_abscissas_.begin(), coarse_abscissas_.end()); | |
512 | for (size_t i = 0; i < samples; ++i) | |
513 | { | |
514 | precise_abscissas_[i] = static_cast<PreciseReal>(coarse_abscissas_[i]); | |
515 | } | |
516 | } | |
517 | ||
518 | precise_ordinates_.resize(samples); | |
519 | for (size_t i = 0; i < samples; ++i) | |
520 | { | |
521 | precise_ordinates_[i] = hi_acc_impl(precise_abscissas_[i]); | |
522 | } | |
523 | ||
524 | cond_.resize(samples, std::numeric_limits<PreciseReal>::quiet_NaN()); | |
525 | for (size_t i = 0 ; i < samples; ++i) | |
526 | { | |
527 | PreciseReal y = precise_ordinates_[i]; | |
528 | if (y != 0) | |
529 | { | |
530 | // Maybe cond_ is badly names; should it be half_cond_? | |
531 | cond_[i] = boost::math::tools::evaluation_condition_number(hi_acc_impl, precise_abscissas_[i])/2; | |
532 | // Half-ULP accuracy is the correctly rounded result, so make sure the envelop doesn't go below this: | |
533 | if (cond_[i] < 0.5) | |
534 | { | |
535 | cond_[i] = 0.5; | |
536 | } | |
537 | } | |
538 | // else leave it as nan. | |
539 | } | |
540 | clip_ = -1; | |
541 | width_ = 1100; | |
542 | envelope_color_ = "chartreuse"; | |
543 | ulp_envelope_ = true; | |
544 | horizontal_lines_ = 8; | |
545 | vertical_lines_ = 10; | |
546 | title_ = ""; | |
547 | background_color_ = "black"; | |
548 | font_color_ = "white"; | |
549 | } | |
550 | ||
551 | template<class F, typename PreciseReal, typename CoarseReal> | |
552 | template<class G> | |
553 | ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::add_fn(G g, std::string const & color) | |
554 | { | |
555 | using std::abs; | |
556 | size_t samples = precise_abscissas_.size(); | |
557 | std::vector<CoarseReal> ulps(samples); | |
558 | for (size_t i = 0; i < samples; ++i) | |
559 | { | |
560 | PreciseReal y_hi_acc = precise_ordinates_[i]; | |
561 | PreciseReal y_lo_acc = static_cast<PreciseReal>(g(coarse_abscissas_[i])); | |
562 | PreciseReal absy = abs(y_hi_acc); | |
563 | PreciseReal dist = static_cast<PreciseReal>(nextafter(static_cast<CoarseReal>(absy), std::numeric_limits<CoarseReal>::max()) - static_cast<CoarseReal>(absy)); | |
564 | ulps[i] = static_cast<CoarseReal>((y_lo_acc - y_hi_acc)/dist); | |
565 | } | |
566 | ulp_list_.emplace_back(ulps); | |
567 | colors_.emplace_back(color); | |
568 | return *this; | |
569 | } | |
570 | ||
571 | ||
572 | ||
573 | ||
574 | } // namespace boost::math::tools | |
575 | #endif |