// (C) Copyright Nick Thompson 2020. // 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) #ifndef BOOST_MATH_TOOLS_ULP_PLOT_HPP #define BOOST_MATH_TOOLS_ULP_PLOT_HPP #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifndef BOOST_MATH_STANDALONE #include #endif // Design of this function comes from: // https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/ // The envelope is the maximum of 1/2 and half the condition number of function evaluation. namespace boost::math::tools { namespace detail { template void write_gridlines(std::ostream& fs, int horizontal_lines, int vertical_lines, F1 x_scale, F2 y_scale, CoarseReal min_x, CoarseReal max_x, PreciseReal min_y, PreciseReal max_y, int graph_width, int graph_height, int margin_left, std::string const & font_color) { // Make a grid: for (int i = 1; i <= horizontal_lines; ++i) { PreciseReal y_cord_dataspace = min_y + ((max_y - min_y)*i)/horizontal_lines; auto y = y_scale(y_cord_dataspace); fs << "\n"; fs << "" << std::setprecision(4) << y_cord_dataspace << "\n"; } for (int i = 1; i <= vertical_lines; ++i) { CoarseReal x_cord_dataspace = min_x + ((max_x - min_x)*i)/vertical_lines; CoarseReal x = x_scale(x_cord_dataspace); fs << "\n"; fs << "" << std::setprecision(4) << x_cord_dataspace << "\n"; } } } template class ulps_plot { public: ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, size_t samples = 1000, bool perturb_abscissas = false, int random_seed = -1); ulps_plot& clip(PreciseReal clip); ulps_plot& width(int width); ulps_plot& envelope_color(std::string const & color); ulps_plot& title(std::string const & title); ulps_plot& background_color(std::string const & background_color); ulps_plot& font_color(std::string const & font_color); ulps_plot& crop_color(std::string const & color); ulps_plot& nan_color(std::string const & color); ulps_plot& ulp_envelope(bool write_ulp); template ulps_plot& add_fn(G g, std::string const & color = "steelblue"); ulps_plot& horizontal_lines(int horizontal_lines); ulps_plot& vertical_lines(int vertical_lines); void write(std::string const & filename) const; friend std::ostream& operator<<(std::ostream& fs, ulps_plot const & plot) { using std::abs; using std::floor; using std::isnan; if (plot.ulp_list_.size() == 0) { throw std::domain_error("No functions added for comparison."); } if (plot.width_ <= 1) { throw std::domain_error("Width = " + std::to_string(plot.width_) + ", which is too small."); } PreciseReal worst_ulp_distance = 0; PreciseReal min_y = (std::numeric_limits::max)(); PreciseReal max_y = std::numeric_limits::lowest(); for (auto const & ulp_vec : plot.ulp_list_) { for (auto const & ulp : ulp_vec) { if (static_cast(abs(ulp)) > worst_ulp_distance) { worst_ulp_distance = static_cast(abs(ulp)); } if (static_cast(ulp) < min_y) { min_y = static_cast(ulp); } if (static_cast(ulp) > max_y) { max_y = static_cast(ulp); } } } // half-ulp accuracy is the best that can be expected; sometimes we can get less, but barely less. // then the axes don't show up; painful! if (max_y < 0.5) { max_y = 0.5; } if (min_y > -0.5) { min_y = -0.5; } if (plot.clip_ > 0) { if (max_y > plot.clip_) { max_y = plot.clip_; } if (min_y < -plot.clip_) { min_y = -plot.clip_; } } int height = static_cast(floor(static_cast(plot.width_)/1.61803)); int margin_top = 40; int margin_left = 25; if (plot.title_.size() == 0) { margin_top = 10; margin_left = 15; } int margin_bottom = 20; int margin_right = 20; int graph_height = height - margin_bottom - margin_top; int graph_width = plot.width_ - margin_left - margin_right; // Maps [a,b] to [0, graph_width] auto x_scale = [&](CoarseReal x)->CoarseReal { return ((x-plot.a_)/(plot.b_ - plot.a_))*static_cast(graph_width); }; auto y_scale = [&](PreciseReal y)->PreciseReal { return ((max_y - y)/(max_y - min_y) )*static_cast(graph_height); }; fs << "\n" << "\n"; return fs; } private: std::vector precise_abscissas_; std::vector coarse_abscissas_; std::vector precise_ordinates_; std::vector cond_; std::list> ulp_list_; std::vector colors_; CoarseReal a_; CoarseReal b_; PreciseReal clip_; int width_; std::string envelope_color_; bool ulp_envelope_; int horizontal_lines_; int vertical_lines_; std::string title_; std::string background_color_; std::string font_color_; std::string crop_color_; std::string nan_color_; }; template ulps_plot& ulps_plot::envelope_color(std::string const & color) { envelope_color_ = color; return *this; } template ulps_plot& ulps_plot::clip(PreciseReal clip) { clip_ = clip; return *this; } template ulps_plot& ulps_plot::width(int width) { width_ = width; return *this; } template ulps_plot& ulps_plot::horizontal_lines(int horizontal_lines) { horizontal_lines_ = horizontal_lines; return *this; } template ulps_plot& ulps_plot::vertical_lines(int vertical_lines) { vertical_lines_ = vertical_lines; return *this; } template ulps_plot& ulps_plot::title(std::string const & title) { title_ = title; return *this; } template ulps_plot& ulps_plot::background_color(std::string const & background_color) { background_color_ = background_color; return *this; } template ulps_plot& ulps_plot::font_color(std::string const & font_color) { font_color_ = font_color; return *this; } template ulps_plot& ulps_plot::crop_color(std::string const & color) { crop_color_ = color; return *this; } template ulps_plot& ulps_plot::nan_color(std::string const & color) { nan_color_ = color; return *this; } template ulps_plot& ulps_plot::ulp_envelope(bool write_ulp_envelope) { ulp_envelope_ = write_ulp_envelope; return *this; } namespace detail{ bool ends_with(std::string const& filename, std::string const& suffix) { if(filename.size() < suffix.size()) { return false; } return std::equal(std::begin(suffix), std::end(suffix), std::end(filename) - suffix.size()); } } template void ulps_plot::write(std::string const & filename) const { if(!boost::math::tools::detail::ends_with(filename, ".svg")) { throw std::logic_error("Only svg files are supported at this time."); } std::ofstream fs(filename); fs << *this; fs.close(); } template ulps_plot::ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, size_t samples, bool perturb_abscissas, int random_seed) : crop_color_("red") { // Use digits10 for this comparison in case the two types have differeing radixes: static_assert(std::numeric_limits::digits10 >= std::numeric_limits::digits10, "PreciseReal must have higher precision that CoarseReal"); if (samples < 10) { throw std::domain_error("Must have at least 10 samples, samples = " + std::to_string(samples)); } if (b <= a) { throw std::domain_error("On interval [a,b], b > a is required."); } a_ = a; b_ = b; std::mt19937_64 gen; if (random_seed == -1) { std::random_device rd; gen.seed(rd()); } // Boost's uniform_real_distribution can generate quad and multiprecision random numbers; std's cannot: #ifndef BOOST_MATH_STANDALONE boost::random::uniform_real_distribution dis(static_cast(a), static_cast(b)); #else // Use std::random in standalone mode if it is a type that the standard library can support (float, double, or long double) static_assert(std::numeric_limits::digits10 <= std::numeric_limits::digits10, "Standalone mode does not support types with precision that exceeds long double"); std::uniform_real_distribution dis(static_cast(a), static_cast(b)); #endif precise_abscissas_.resize(samples); coarse_abscissas_.resize(samples); if (perturb_abscissas) { for(size_t i = 0; i < samples; ++i) { precise_abscissas_[i] = dis(gen); } std::sort(precise_abscissas_.begin(), precise_abscissas_.end()); for (size_t i = 0; i < samples; ++i) { coarse_abscissas_[i] = static_cast(precise_abscissas_[i]); } } else { for(size_t i = 0; i < samples; ++i) { coarse_abscissas_[i] = static_cast(dis(gen)); } std::sort(coarse_abscissas_.begin(), coarse_abscissas_.end()); for (size_t i = 0; i < samples; ++i) { precise_abscissas_[i] = static_cast(coarse_abscissas_[i]); } } precise_ordinates_.resize(samples); for (size_t i = 0; i < samples; ++i) { precise_ordinates_[i] = hi_acc_impl(precise_abscissas_[i]); } cond_.resize(samples, std::numeric_limits::quiet_NaN()); for (size_t i = 0 ; i < samples; ++i) { PreciseReal y = precise_ordinates_[i]; if (y != 0) { // Maybe cond_ is badly names; should it be half_cond_? cond_[i] = boost::math::tools::evaluation_condition_number(hi_acc_impl, precise_abscissas_[i])/2; // Half-ULP accuracy is the correctly rounded result, so make sure the envelop doesn't go below this: if (cond_[i] < 0.5) { cond_[i] = 0.5; } } // else leave it as nan. } clip_ = -1; width_ = 1100; envelope_color_ = "chartreuse"; ulp_envelope_ = true; horizontal_lines_ = 8; vertical_lines_ = 10; title_ = ""; background_color_ = "black"; font_color_ = "white"; } template template ulps_plot& ulps_plot::add_fn(G g, std::string const & color) { using std::abs; size_t samples = precise_abscissas_.size(); std::vector ulps(samples); for (size_t i = 0; i < samples; ++i) { PreciseReal y_hi_acc = precise_ordinates_[i]; PreciseReal y_lo_acc = static_cast(g(coarse_abscissas_[i])); PreciseReal absy = abs(y_hi_acc); PreciseReal dist = static_cast(nextafter(static_cast(absy), (std::numeric_limits::max)()) - static_cast(absy)); ulps[i] = static_cast((y_lo_acc - y_hi_acc)/dist); } ulp_list_.emplace_back(ulps); colors_.emplace_back(color); return *this; } } // namespace boost::math::tools #endif