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