ceph/src/boost/libs/math/example/daubechies_wavelets/daubechies_wavelet_plots.cpp

   1 /*
   2  * Copyright Nick Thompson, 2020
   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 #include <iostream>
   8 #include <boost/core/demangle.hpp>
   9 #include <boost/hana/for_each.hpp>
  10 #include <boost/hana/ext/std/integer_sequence.hpp>
  11
  12 #include <boost/multiprecision/float128.hpp>
  13 #include <boost/math/special_functions/daubechies_wavelet.hpp>
  14 #include <quicksvg/graph_fn.hpp>
  15 #include <quicksvg/ulp_plot.hpp>
  16
  17
  18 using boost::multiprecision::float128;
  19 constexpr const int GRAPH_WIDTH = 700;
  20
  21 template<typename Real, int p>
  22 void plot_psi(int grid_refinements = -1)
  23 {
  24     auto psi = boost::math::daubechies_wavelet<Real, p>();
  25     if (grid_refinements >= 0)
  26     {
  27         psi = boost::math::daubechies_wavelet<Real, p>(grid_refinements);
  28     }
  29     auto [a, b] = psi.support();
  30     std::string title = "Daubechies " + std::to_string(p) + " wavelet";
  31     title = "";
  32     std::string filename = "daubechies_" + std::to_string(p) + "_wavelet.svg";
  33     int samples = 1024;
  34     quicksvg::graph_fn daub(a, b, title, filename, samples, GRAPH_WIDTH);
  35     daub.set_gridlines(8, 2*p-1);
  36     daub.set_stroke_width(1);
  37     daub.add_fn(psi);
  38     daub.write_all();
  39 }
  40
  41 template<typename Real, int p>
  42 void plot_dpsi(int grid_refinements = -1)
  43 {
  44     auto psi = boost::math::daubechies_wavelet<Real, p>();
  45     if (grid_refinements >= 0)
  46     {
  47         psi = boost::math::daubechies_wavelet<Real, p>(grid_refinements);
  48     }
  49     auto [a, b] = psi.support();
  50     std::string title = "Daubechies " + std::to_string(p) + " wavelet derivative";
  51     title = "";
  52     std::string filename = "daubechies_" + std::to_string(p) + "_wavelet_prime.svg";
  53     int samples = 1024;
  54     quicksvg::graph_fn daub(a, b, title, filename, samples, GRAPH_WIDTH);
  55     daub.set_stroke_width(1);
  56     daub.set_gridlines(8, 2*p-1);
  57     auto dpsi = [psi](Real x)->Real { return psi.prime(x); };
  58     daub.add_fn(dpsi);
  59     daub.write_all();
  60 }
  61
  62 template<typename Real, int p>
  63 void plot_convergence()
  64 {
  65     auto psi1 = boost::math::daubechies_wavelet<Real, p>(1);
  66     auto [a, b] = psi1.support();
  67     std::string title = "Daubechies " + std::to_string(p) + " wavelet at  1 (orange), 2 (red), and 21 (blue) grid refinements";
  68     title = "";
  69     std::string filename = "daubechies_" + std::to_string(p) + "_wavelet_convergence.svg";
  70
  71     quicksvg::graph_fn daub(a, b, title, filename, 1024, GRAPH_WIDTH);
  72     daub.set_stroke_width(1);
  73     daub.set_gridlines(8, 2*p-1);
  74
  75     daub.add_fn(psi1, "orange");
  76     auto psi2 = boost::math::daubechies_wavelet<Real, p>(2);
  77     daub.add_fn(psi2, "red");
  78
  79     auto psi21 = boost::math::daubechies_wavelet<Real, p>(21);
  80     daub.add_fn(psi21);
  81
  82     daub.write_all();
  83 }
  84
  85 template<typename Real, int p>
  86 void plot_condition_number()
  87 {
  88     using std::abs;
  89     using std::log;
  90     static_assert(p >= 3, "p = 2 is not differentiable, so condition numbers cannot be effectively evaluated.");
  91     auto phi = boost::math::daubechies_wavelet<Real, p>();
  92     Real a = phi.support().first + 1000*std::sqrt(std::numeric_limits<Real>::epsilon());
  93     Real b = phi.support().second - 1000*std::sqrt(std::numeric_limits<Real>::epsilon());
  94     std::string title = "log10 of condition number of function evaluation for Daubechies " + std::to_string(p) + " wavelet function.";
  95     title = "";
  96     std::string filename = "daubechies_" + std::to_string(p) + "_wavelet_condition_number.svg";
  97
  98
  99     quicksvg::graph_fn daub(a, b, title, filename, 2048, GRAPH_WIDTH);
 100     daub.set_stroke_width(1);
 101     daub.set_gridlines(8, 2*p-1);
 102
 103     auto cond = [&phi](Real x)
 104     {
 105         Real y = phi(x);
 106         Real dydx = phi.prime(x);
 107         Real z = abs(x*dydx/y);
 108         using std::isnan;
 109         if (z==0)
 110         {
 111             return Real(-1);
 112         }
 113         if (isnan(z))
 114         {
 115             // Graphing libraries don't like nan's:
 116             return Real(1);
 117         }
 118         return log10(z);
 119     };
 120     daub.add_fn(cond);
 121     daub.write_all();
 122 }
 123
 124 template<typename CoarseReal, typename PreciseReal, int p, class PsiPrecise>
 125 void do_ulp(int coarse_refinements, PsiPrecise psi_precise)
 126 {
 127     auto psi_coarse = boost::math::daubechies_wavelet<CoarseReal, p>(coarse_refinements);
 128
 129     std::string title = std::to_string(p) + " vanishing moment ULP plot at " + std::to_string(coarse_refinements) + " refinements and " + boost::core::demangle(typeid(CoarseReal).name()) + " precision";
 130     title = "";
 131
 132     std::string filename = "daubechies_" + std::to_string(p) + "_wavelet_" + boost::core::demangle(typeid(CoarseReal).name()) + "_" + std::to_string(coarse_refinements) + "_refinements.svg";
 133     int samples = 20000;
 134     int clip = 20;
 135     int horizontal_lines = 8;
 136     int vertical_lines = 2*p - 1;
 137     quicksvg::ulp_plot<decltype(psi_coarse), CoarseReal, decltype(psi_precise), PreciseReal>(psi_coarse, psi_precise, CoarseReal(psi_coarse.support().first), psi_coarse.support().second, title, filename, samples, GRAPH_WIDTH, clip, horizontal_lines, vertical_lines);
 138 }
 139
 140
 141 int main()
 142 {
 143     boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i){ plot_psi<double, i+2>(); });
 144     boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i){ plot_dpsi<double, i+3>(); });
 145     boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i){ plot_condition_number<double, i+3>(); });
 146     boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i){ plot_convergence<double, i+2>(); });
 147
 148     using PreciseReal = float128;
 149     using CoarseReal = double;
 150     int precise_refinements = 22;
 151     constexpr const int p = 9;
 152     std::cout << "Computing precise wavelet function in " << boost::core::demangle(typeid(PreciseReal).name()) << " precision.\n";
 153     auto phi_precise = boost::math::daubechies_wavelet<PreciseReal, p>(precise_refinements);
 154     std::cout << "Beginning comparison with functions computed in " << boost::core::demangle(typeid(CoarseReal).name()) << " precision.\n";
 155     for (int i = 7; i <= precise_refinements-1; ++i)
 156     {
 157         std::cout << "\tCoarse refinement " << i << "\n";
 158         do_ulp<CoarseReal, PreciseReal, p>(i, phi_precise);
 159     }
 160 }