1 /*
2 * Copyright Nick Thompson, John Maddock 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
8 #define BOOST_MATH_GENERATE_DAUBECHIES_GRID
9
10 #include <iostream>
11 #include <vector>
12 #include <numeric>
13 #include <list>
14 #include <cmath>
15 #include <cassert>
16 #include <fstream>
17 #include <Eigen/Eigenvalues>
18 #include <boost/hana/for_each.hpp>
19 #include <boost/hana/ext/std/integer_sequence.hpp>
20 #include <boost/core/demangle.hpp>
21 #ifdef BOOST_HAS_FLOAT128
22 #include <boost/multiprecision/float128.hpp>
23 #endif
24 #include <boost/math/constants/constants.hpp>
25 #include <boost/math/filters/daubechies.hpp>
26 #include <boost/math/special_functions/factorials.hpp>
27 #include <boost/multiprecision/cpp_bin_float.hpp>
28
29 typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<237, boost::multiprecision::backends::digit_base_2, std::allocator<char>, boost::int32_t, -262142, 262143>, boost::multiprecision::et_off> octuple_type;
30
31 #ifdef BOOST_HAS_FLOAT128
32 typedef boost::multiprecision::float128 float128_t;
33 #else
34 typedef boost::multiprecision::cpp_bin_float_quad float128_t;
35 #endif
36
37 template<class Real, int p>
integer_grid()38 std::list<std::vector<Real>> integer_grid()
39 {
40 std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 3);
41 using std::abs;
42 using std::sqrt;
43 using std::pow;
44 std::list<std::vector<Real>> grids;
45
46 auto c = boost::math::filters::daubechies_scaling_filter<Real, p>();
47 for (auto & x : c)
48 {
49 x *= boost::math::constants::root_two<Real>();
50 }
51 std::cout << "\n\nTaps in filter = " << c.size() << "\n";
52
53
54 Eigen::Matrix<Real, 2*p - 2, 2*p-2> A;
55 for (int j = 0; j < 2*p-2; ++j) {
56 for (int k = 0; k < 2*p-2; ++k) {
57 if ( (2*j-k + 1) < 0 || (2*j - k + 1) >= 2*p)
58 {
59 A(j,k) = 0;
60 }
61 else {
62 A(j,k) = c[2*j - k + 1];
63 }
64 }
65 }
66
67 Eigen::EigenSolver<decltype(A)> es(A);
68
69 auto complex_eigs = es.eigenvalues();
70
71 std::vector<Real> eigs(complex_eigs.size(), std::numeric_limits<Real>::quiet_NaN());
72
73 std::cout << "Eigenvalues = {";
74 for (long i = 0; i < complex_eigs.size(); ++i) {
75 assert(abs(complex_eigs[i].imag()) < std::numeric_limits<Real>::epsilon());
76 eigs[i] = complex_eigs[i].real();
77 std::cout << eigs[i] << ", ";
78 }
79 std::cout << "}\n";
80
81 // Eigen does not sort the eigenpairs by any criteria on the eigenvalues.
82 // In any case, even if it did, some of the eigenpairs do not correspond to derivatives anyway.
83 for (size_t j = 0; j < eigs.size(); ++j) {
84 auto f = [&](Real x) {
85 return abs(x - Real(1)/Real(1 << j) ) < sqrt(std::numeric_limits<Real>::epsilon());
86 };
87 auto it = std::find_if(eigs.begin(), eigs.end(), f);
88 if (it == eigs.end()) {
89 std::cout << "couldn't find eigenvalue " << Real(1)/Real(1 << j) << "\n";
90 continue;
91 }
92 size_t idx = std::distance(eigs.begin(), it);
93 std::cout << "Eigenvector for derivative " << j << " is at index " << idx << "\n";
94 auto eigenvector_matrix = es.eigenvectors();
95 auto complex_eigenvec = eigenvector_matrix.col(idx);
96
97 std::vector<Real> eigenvec(complex_eigenvec.size() + 2, std::numeric_limits<Real>::quiet_NaN());
98 eigenvec[0] = 0;
99 eigenvec[eigenvec.size()-1] = 0;
100 for (size_t i = 0; i < eigenvec.size() - 2; ++i) {
101 assert(abs(complex_eigenvec[i].imag()) < std::numeric_limits<Real>::epsilon());
102 eigenvec[i+1] = complex_eigenvec[i].real();
103 }
104
105 Real sum = 0;
106 for(size_t k = 1; k < eigenvec.size(); ++k) {
107 sum += pow(k, j)*eigenvec[k];
108 }
109
110 Real alpha = pow(-1, j)*boost::math::factorial<Real>(j)/sum;
111
112 for (size_t i = 1; i < eigenvec.size(); ++i) {
113 eigenvec[i] *= alpha;
114 }
115
116
117 std::cout << "Eigenvector = {";
118 for (size_t i = 0; i < eigenvec.size() -1; ++i) {
119 std::cout << eigenvec[i] << ", ";
120 }
121 std::cout << eigenvec[eigenvec.size()-1] << "}\n";
122
123 sum = 0;
124 for(size_t k = 1; k < eigenvec.size(); ++k) {
125 sum += pow(k, j)*eigenvec[k];
126 }
127
128 std::cout << "Moment sum = " << sum << ", expected = " << pow(-1, j)*boost::math::factorial<Real>(j) << "\n";
129
130 assert(abs(sum - pow(-1, j)*boost::math::factorial<Real>(j))/abs(pow(-1, j)*boost::math::factorial<Real>(j)) < sqrt(std::numeric_limits<Real>::epsilon()));
131
132 grids.push_back(eigenvec);
133 }
134
135
136 return grids;
137 }
138
139 template<class Real, int p>
write_grid(std::ofstream & fs)140 void write_grid(std::ofstream & fs)
141 {
142 auto grids = integer_grid<Real, p>();
143 size_t j = 0;
144 fs << std::setprecision(std::numeric_limits< boost::multiprecision::cpp_bin_float_quad>::max_digits10);
145 for (auto it = grids.begin(); it != grids.end(); ++it)
146 {
147 auto const& grid = *it;
148 fs << "template <typename Real> struct daubechies_scaling_integer_grid_imp <Real, " << p << ", ";
149 fs << j << "> { static inline constexpr std::array<Real, " << grid.size() << "> value = { ";
150 for (size_t i = 0; i < grid.size() -1; ++i){
151 fs << "C_(" << static_cast<float128_t>(grid[i]) << "), ";
152 }
153 fs << "C_(" << static_cast<float128_t>(grid[grid.size()-1]) << ") }; };\n";
154 ++j;
155 }
156 }
157
main()158 int main()
159 {
160 constexpr const size_t p_max = 18;
161 std::ofstream fs{"daubechies_scaling_integer_grid.hpp"};
162 fs << "/*\n"
163 << " * Copyright Nick Thompson, John Maddock 2020\n"
164 << " * Use, modification and distribution are subject to the\n"
165 << " * Boost Software License, Version 1.0. (See accompanying file\n"
166 << " * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)\n"
167 << " */\n"
168 << "// THIS FILE GENERATED BY EXAMPLE/DAUBECHIES_SCALING_INTEGER_GRID.CPP, DO NOT EDIT.\n"
169 << "#ifndef BOOST_MATH_DAUBECHIES_SCALING_INTEGER_GRID_HPP\n"
170 << "#define BOOST_MATH_DAUBECHIES_SCALING_INTEGER_GRID_HPP\n"
171 << "#include <array>\n"
172 << "#include <float.h>\n"
173 << "#include <boost/config.hpp>\n"
174 << "/*\n"
175 << "In order to keep the character count as small as possible and speed up\n"
176 << "compiler parsing times, we define a macro C_ which appends an appropriate\n"
177 << "suffix to each literal, and then casts it to type Real.\n"
178 << "The suffix is as follows:\n\n"
179 << "* Q, when we have __float128 support.\n"
180 << "* L, when we have either 80 or 128 bit long doubles.\n"
181 << "* Nothing otherwise.\n"
182 << "*/\n\n"
183 << "#ifdef BOOST_HAS_FLOAT128\n"
184 << "# define C_(x) static_cast<Real>(x##Q)\n"
185 << "#elif (LDBL_MANT_DIG > DBL_MANT_DIG)\n"
186 << "# define C_(x) static_cast<Real>(x##L)\n"
187 << "#else\n"
188 << "# define C_(x) static_cast<Real>(x)\n"
189 << "#endif\n\n"
190 << "namespace boost::math::detail {\n\n"
191 << "template <typename Real, int p, int order> struct daubechies_scaling_integer_grid_imp;\n\n";
192
193 fs << std::hexfloat << std::setprecision(std::numeric_limits<boost::multiprecision::cpp_bin_float_quad>::max_digits10);
194
195 boost::hana::for_each(std::make_index_sequence<p_max>(), [&](auto idx){
196 write_grid<octuple_type, idx+2>(fs);
197 });
198
199 fs << "\n\ntemplate <typename Real, unsigned p, unsigned order>\n"
200 << "constexpr inline std::array<Real, 2*p> daubechies_scaling_integer_grid()\n"
201 << "{\n"
202 << " static_assert(sizeof(Real) <= 16, \"Integer grids only computed up to 128 bits of precision.\");\n"
203 << " static_assert(p <= " << p_max + 1 << ", \"Integer grids only implemented up to " << p_max + 1 << ".\");\n"
204 << " static_assert(p > 1, \"Integer grids only implemented for p >= 2.\");\n"
205 << " return daubechies_scaling_integer_grid_imp<Real, p, order>::value;\n"
206 << "}\n\n";
207
208 fs << "} // namespaces\n";
209 fs << "#endif\n";
210 fs.close();
211
212 return 0;
213 }
214