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>
plot_psi(int grid_refinements=-1)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>
plot_dpsi(int grid_refinements=-1)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>
plot_convergence()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>
plot_condition_number()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>
do_ulp(int coarse_refinements,PsiPrecise psi_precise)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
main()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 }
161