Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : The VES code module is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "WaveletGrid.h"
24 : #include "tools/Exception.h"
25 : #include "tools/Grid.h"
26 : #include "tools/Matrix.h"
27 :
28 : #include <memory>
29 : #include <unordered_map>
30 : #include <vector>
31 :
32 : namespace PLMD {
33 : namespace ves {
34 :
35 :
36 : // construction of Wavelet grid according to the Daubechies-Lagarias method
37 : // see Strang, Nguyen "Wavelets and Filter Banks" chapter 6.3
38 47 : std::unique_ptr<Grid> WaveletGrid::setupGrid(const unsigned order, unsigned gridsize, const bool use_mother_wavelet, const Type type) {
39 47 : plumed_assert(order>=1) << "Wavelet order has to be a positive integer";
40 : // calculate the grid properties of the scaling grid
41 : // the range of the grid is from 0 to maxsupport
42 47 : unsigned maxsupport = order*2 -1;
43 : // determine needed recursion depth for specified size
44 : unsigned recursion_number = 0;
45 173 : while (maxsupport*(1U<<recursion_number) < gridsize) {
46 126 : recursion_number++;
47 : }
48 : // set "true" gridsize
49 47 : unsigned bins_per_int = 1U<<recursion_number;
50 : gridsize = maxsupport*bins_per_int;
51 :
52 : // Filter coefficients
53 47 : std::vector<double> h_coeffs = getFilterCoefficients(order, true, type);
54 : // Vector with the Matrices M0 and M1 for the cascade
55 47 : std::vector<Matrix<double>> h_Matvec = setupMatrices(h_coeffs);
56 : std::vector<Matrix<double>> g_Matvec; // only filled if needed for wavelet
57 :
58 : // get the values at integers
59 47 : std::vector<double> values_at_integers = calcIntegerValues(h_Matvec[0], 0);
60 47 : std::vector<double> derivs_at_integers = calcIntegerValues(h_Matvec[0], 1);
61 :
62 : std::string gridvarname; // stores the name of the grid variable
63 47 : if (use_mother_wavelet) { // get the highpass filter coefficients as well
64 2 : std::vector<double> g_coeffs = getFilterCoefficients(order, false, type);
65 4 : g_Matvec = setupMatrices(g_coeffs);
66 4 : gridvarname = typeToString(type,true)+std::to_string(order)+"_psi";
67 : } else {
68 90 : gridvarname = typeToString(type,true)+std::to_string(order)+"_phi";
69 : }
70 :
71 : // Set up the grid with correct properties
72 : auto grid = std::unique_ptr<Grid>(new Grid(gridvarname, {"position"}, {"0"},
73 517 : {std::to_string(maxsupport)}, {gridsize}, false, true, {false}, {"0."}, {"0."}));
74 :
75 47 : BinaryMap values = cascade(h_Matvec, g_Matvec, values_at_integers, recursion_number, bins_per_int, 0, use_mother_wavelet);
76 47 : BinaryMap derivs = cascade(h_Matvec, g_Matvec, derivs_at_integers, recursion_number, bins_per_int, 1, use_mother_wavelet);
77 :
78 47 : fillGridFromMaps(grid, values, derivs);
79 :
80 47 : return grid;
81 47 : }
82 :
83 :
84 :
85 49 : std::vector<Matrix<double>> WaveletGrid::setupMatrices(const std::vector<double>& coeffs) {
86 : Matrix<double> M0, M1;
87 49 : const int N = coeffs.size() -1;
88 : M0.resize(N,N);
89 : M1.resize(N,N);
90 1108 : for (int i = 0; i < N; ++i) {
91 27644 : for (int j = 0; j < N; ++j) {
92 26585 : int shift = 2*i -j;
93 26585 : if (0 <= shift && shift <= N) {
94 13822 : M0[i][j] = coeffs[2*i -j]; // normalization is already included in coeffs
95 : }
96 26585 : if (-1 <= shift && shift <= N -1) {
97 13822 : M1[i][j] = coeffs[2*i -j + 1];
98 : }
99 : }
100 : }
101 196 : return std::vector<Matrix<double>> {M0, M1};
102 : }
103 :
104 :
105 94 : std::vector<double> WaveletGrid::calcIntegerValues(const Matrix<double> &M, const int deriv) {
106 : // corresponding eigenvalue of the matrix
107 : double eigenvalue = std::pow(0.5, deriv);
108 94 : std::vector<double> values = getEigenvector(M, eigenvalue);
109 :
110 : // normalization of the eigenvector
111 : double normfactor = 0.;
112 : // i=0 is always 0; for derivs > 1 an additional factorial would have to be added
113 2042 : for (unsigned i=1; i<values.size(); ++i) {
114 1948 : normfactor += values[i] * std::pow(-i, deriv);
115 : }
116 94 : normfactor = 1/normfactor;
117 2136 : for (auto& value : values) {
118 2042 : value *= normfactor;
119 : }
120 :
121 94 : return values;
122 : }
123 :
124 :
125 : // maybe move this to the tools/matrix.h file?
126 : // this works reliably only for singular eigenvalues
127 : //
128 94 : std::vector<double> WaveletGrid::getEigenvector(const Matrix<double> &A, const double eigenvalue) {
129 : // mostly copied from tools/matrix.h
130 94 : int info, N = A.ncols(); // ncols == nrows
131 94 : std::vector<double> da(N*N);
132 94 : std::vector<double> S(N);
133 94 : std::vector<double> U(N*N);
134 94 : std::vector<double> VT(N*N);
135 94 : std::vector<int> iwork(8*N);
136 :
137 : // Transfer the matrix to the local array and substract eigenvalue
138 2136 : for (int i=0; i<N; ++i)
139 53768 : for (int j=0; j<N; ++j) {
140 51726 : da[i*N+j]=static_cast<double>( A(j,i) );
141 51726 : if (i==j) {
142 2042 : da[i*N+j] -= eigenvalue;
143 : }
144 : }
145 :
146 : // This optimizes the size of the work array used in lapack singular value decomposition
147 94 : int lwork=-1;
148 94 : std::vector<double> work(1);
149 94 : plumed_lapack_dgesdd( "A", &N, &N, da.data(), &N, S.data(), U.data(), &N, VT.data(), &N, work.data(), &lwork, iwork.data(), &info );
150 :
151 : // Retrieve correct sizes for work and reallocate
152 94 : lwork=static_cast<int>(work[0]);
153 94 : work.resize(lwork);
154 :
155 : // This does the singular value decomposition
156 94 : plumed_lapack_dgesdd( "A", &N, &N, da.data(), &N, S.data(), U.data(), &N, VT.data(), &N, work.data(), &lwork, iwork.data(), &info );
157 :
158 : // fill eigenvector with last column of VT
159 94 : std::vector<double> eigenvector(N);
160 2136 : for (int i=0; i<N; ++i) {
161 2042 : eigenvector[i] = VT[N-1 + i*N];
162 : }
163 :
164 94 : return eigenvector;
165 : }
166 :
167 :
168 94 : WaveletGrid::BinaryMap WaveletGrid::cascade(std::vector<Matrix<double>>& h_Matvec, std::vector<Matrix<double>>& g_Matvec, const std::vector<double>& values_at_integers, const unsigned recursion_number, const unsigned bins_per_int, const unsigned derivnum, const bool use_mother_wavelet) {
169 : BinaryMap scaling_map, wavelet_map;
170 94 : scaling_map.reserve(bins_per_int);
171 : // vector to store the binary representation of all the decimal parts
172 : std::vector<std::string> binary_vec;
173 : // vector used as result of the matrix multiplications
174 : std::vector<double> temp_values;
175 :
176 : // multiply matrices by 2 if derivatives are calculated (assumes ascending order)
177 94 : if (derivnum != 0)
178 141 : for (auto& M : h_Matvec) {
179 188 : M *= 2;
180 : }
181 :
182 94 : if (use_mother_wavelet) {
183 : wavelet_map.reserve(bins_per_int);
184 4 : if (derivnum != 0)
185 6 : for (auto& M : g_Matvec) {
186 8 : M *= 2;
187 : }
188 : }
189 :
190 : // fill the first two datasets by hand
191 282 : scaling_map["0"] = values_at_integers;
192 94 : mult(h_Matvec[1], values_at_integers, temp_values);
193 188 : scaling_map["1"] = temp_values;
194 :
195 94 : if (use_mother_wavelet) {
196 4 : mult(g_Matvec[0], values_at_integers, temp_values);
197 12 : wavelet_map["0"] = temp_values;
198 4 : mult(g_Matvec[1], values_at_integers, temp_values);
199 12 : wavelet_map["1"] = temp_values;
200 : }
201 :
202 : // now do the cascade
203 94 : binary_vec.emplace_back("1");
204 252 : for (unsigned i=1; i<recursion_number; ++i) {
205 : std::vector<std::string> new_binary_vec;
206 1144 : for (const auto& binary : binary_vec) {
207 2958 : for (unsigned k=0; k<2; ++k) {
208 : // prepend the new bit
209 3944 : std::string new_binary = std::to_string(k) + binary;
210 1972 : mult(h_Matvec[k], scaling_map[binary], temp_values);
211 1972 : scaling_map.insert(std::pair<std::string, std::vector<double>>(new_binary, temp_values));
212 1972 : if (use_mother_wavelet) {
213 8 : mult(g_Matvec[k], scaling_map[binary], temp_values);
214 16 : wavelet_map.insert(std::pair<std::string, std::vector<double>>(new_binary, temp_values));
215 : }
216 1972 : new_binary_vec.push_back(new_binary);
217 : }
218 : }
219 158 : binary_vec = new_binary_vec;
220 158 : }
221 :
222 188 : return use_mother_wavelet ? wavelet_map : scaling_map;
223 94 : }
224 :
225 :
226 : // Fill the Grid with the values of the unordered maps
227 47 : void WaveletGrid::fillGridFromMaps(std::unique_ptr<Grid>& grid, const BinaryMap& values_map, const BinaryMap& derivs_map) {
228 47 : unsigned bins_per_int = values_map.size();
229 : // this is somewhat complicated… not sure if the unordered_map way is the best way for c++
230 1127 : for (const auto& value_iter : values_map) {
231 : // get decimal of binary key
232 1080 : unsigned decimal = std::stoi(value_iter.first, nullptr, 2);
233 : // corresponding iterator of deriv
234 1080 : auto deriv_iter = derivs_map.find(value_iter.first);
235 : // calculate first grid element
236 1080 : unsigned first_grid_element = decimal * (bins_per_int >> value_iter.first.length());
237 19480 : for (unsigned i=0; i<value_iter.second.size(); ++i) {
238 : // derivative has to be passed as vector
239 18400 : std::vector<double> deriv {deriv_iter->second[i]};
240 18400 : grid->setValueAndDerivatives(first_grid_element + bins_per_int*i, value_iter.second[i], deriv);
241 : }
242 : }
243 47 : }
244 :
245 :
246 47 : WaveletGrid::Type WaveletGrid::stringToType(std::string& type_str) {
247 188 : std::unordered_map<std::string, Type> typemap = { {"DAUBECHIES", Type::db}, {"SYMLETS", Type::sym} };
248 : try {
249 47 : return typemap.at(type_str);
250 0 : } catch(const std::out_of_range& e) {
251 0 : plumed_merror("The specified wavelet type "+type_str+" is not implemented.");
252 0 : }
253 : }
254 :
255 :
256 47 : std::string WaveletGrid::typeToString(Type type, bool abbrev) {
257 : std::string type_str;
258 47 : switch(type) {
259 23 : case Type::db:
260 23 : type_str = abbrev ? "Db" : "DAUBECHIES";
261 : break;
262 24 : case Type::sym:
263 24 : type_str = abbrev ? "Sym" : "SYMLETS";
264 : break;
265 : }
266 47 : return type_str;
267 : }
268 :
269 :
270 : }
271 : }
|
