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