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) { recursion_number++; }
46 : // set "true" gridsize
47 47 : unsigned bins_per_int = 1U<<recursion_number;
48 : gridsize = maxsupport*bins_per_int;
49 :
50 : // Filter coefficients
51 47 : std::vector<double> h_coeffs = getFilterCoefficients(order, true, type);
52 : // Vector with the Matrices M0 and M1 for the cascade
53 47 : std::vector<Matrix<double>> h_Matvec = setupMatrices(h_coeffs);
54 : std::vector<Matrix<double>> g_Matvec; // only filled if needed for wavelet
55 :
56 : // get the values at integers
57 47 : std::vector<double> values_at_integers = calcIntegerValues(h_Matvec[0], 0);
58 47 : std::vector<double> derivs_at_integers = calcIntegerValues(h_Matvec[0], 1);
59 :
60 : std::string gridvarname; // stores the name of the grid variable
61 47 : if (use_mother_wavelet) { // get the highpass filter coefficients as well
62 2 : std::vector<double> g_coeffs = getFilterCoefficients(order, false, type);
63 4 : g_Matvec = setupMatrices(g_coeffs);
64 4 : gridvarname = typeToString(type,true)+std::to_string(order)+"_psi";
65 : }
66 : else {
67 90 : gridvarname = typeToString(type,true)+std::to_string(order)+"_phi";
68 : }
69 :
70 : // Set up the grid with correct properties
71 : auto grid = std::unique_ptr<Grid>(new Grid(gridvarname, {"position"}, {"0"},
72 564 : {std::to_string(maxsupport)}, {gridsize}, false, true, {false}, {"0."}, {"0."}));
73 :
74 47 : BinaryMap values = cascade(h_Matvec, g_Matvec, values_at_integers, recursion_number, bins_per_int, 0, use_mother_wavelet);
75 47 : BinaryMap derivs = cascade(h_Matvec, g_Matvec, derivs_at_integers, recursion_number, bins_per_int, 1, use_mother_wavelet);
76 :
77 47 : fillGridFromMaps(grid, values, derivs);
78 :
79 47 : return grid;
80 47 : }
81 :
82 :
83 :
84 49 : std::vector<Matrix<double>> WaveletGrid::setupMatrices(const std::vector<double>& coeffs) {
85 : Matrix<double> M0, M1;
86 49 : const int N = coeffs.size() -1;
87 : M0.resize(N,N); M1.resize(N,N);
88 1108 : for (int i = 0; i < N; ++i) {
89 27644 : for (int j = 0; j < N; ++j) {
90 26585 : int shift = 2*i -j;
91 26585 : if (0 <= shift && shift <= N) {
92 13822 : M0[i][j] = coeffs[2*i -j]; // normalization is already included in coeffs
93 : }
94 26585 : if (-1 <= shift && shift <= N -1) {
95 13822 : M1[i][j] = coeffs[2*i -j + 1];
96 : }
97 : }
98 : }
99 196 : return std::vector<Matrix<double>> {M0, M1};
100 : }
101 :
102 :
103 94 : std::vector<double> WaveletGrid::calcIntegerValues(const Matrix<double> &M, const int deriv) {
104 : // corresponding eigenvalue of the matrix
105 : double eigenvalue = std::pow(0.5, deriv);
106 94 : std::vector<double> values = getEigenvector(M, eigenvalue);
107 :
108 : // normalization of the eigenvector
109 : double normfactor = 0.;
110 : // i=0 is always 0; for derivs > 1 an additional factorial would have to be added
111 2042 : for (unsigned i=1; i<values.size(); ++i) {
112 1948 : normfactor += values[i] * std::pow(-i, deriv);
113 : }
114 94 : normfactor = 1/normfactor;
115 2136 : for (auto& value : values) {
116 2042 : value *= normfactor;
117 : }
118 :
119 94 : return values;
120 : }
121 :
122 :
123 : // maybe move this to the tools/matrix.h file?
124 : // this works reliably only for singular eigenvalues
125 : //
126 94 : std::vector<double> WaveletGrid::getEigenvector(const Matrix<double> &A, const double eigenvalue) {
127 : // mostly copied from tools/matrix.h
128 94 : int info, N = A.ncols(); // ncols == nrows
129 94 : std::vector<double> da(N*N);
130 94 : std::vector<double> S(N);
131 94 : std::vector<double> U(N*N);
132 94 : std::vector<double> VT(N*N);
133 94 : std::vector<int> iwork(8*N);
134 :
135 : // Transfer the matrix to the local array and substract eigenvalue
136 53862 : for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) {
137 51726 : da[i*N+j]=static_cast<double>( A(j,i) );
138 51726 : if (i==j) da[i*N+j] -= eigenvalue;
139 : }
140 :
141 : // This optimizes the size of the work array used in lapack singular value decomposition
142 94 : int lwork=-1;
143 94 : std::vector<double> work(1);
144 94 : plumed_lapack_dgesdd( "A", &N, &N, da.data(), &N, S.data(), U.data(), &N, VT.data(), &N, work.data(), &lwork, iwork.data(), &info );
145 :
146 : // Retrieve correct sizes for work and reallocate
147 94 : lwork=static_cast<int>(work[0]); work.resize(lwork);
148 :
149 : // This does the singular value decomposition
150 94 : plumed_lapack_dgesdd( "A", &N, &N, da.data(), &N, S.data(), U.data(), &N, VT.data(), &N, work.data(), &lwork, iwork.data(), &info );
151 :
152 : // fill eigenvector with last column of VT
153 94 : std::vector<double> eigenvector(N);
154 2136 : for (int i=0; i<N; ++i) eigenvector[i] = VT[N-1 + i*N];
155 :
156 94 : return eigenvector;
157 : }
158 :
159 :
160 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) {
161 : BinaryMap scaling_map, wavelet_map;
162 94 : scaling_map.reserve(bins_per_int);
163 : // vector to store the binary representation of all the decimal parts
164 : std::vector<std::string> binary_vec;
165 : // vector used as result of the matrix multiplications
166 : std::vector<double> temp_values;
167 :
168 : // multiply matrices by 2 if derivatives are calculated (assumes ascending order)
169 188 : if (derivnum != 0) for (auto& M : h_Matvec) M *= 2;
170 :
171 94 : if (use_mother_wavelet) {
172 : wavelet_map.reserve(bins_per_int);
173 8 : if (derivnum != 0) for (auto& M : g_Matvec) M *= 2;
174 : }
175 :
176 : // fill the first two datasets by hand
177 282 : scaling_map["0"] = values_at_integers;
178 94 : mult(h_Matvec[1], values_at_integers, temp_values);
179 188 : scaling_map["1"] = temp_values;
180 :
181 94 : if (use_mother_wavelet) {
182 4 : mult(g_Matvec[0], values_at_integers, temp_values);
183 12 : wavelet_map["0"] = temp_values;
184 4 : mult(g_Matvec[1], values_at_integers, temp_values);
185 12 : wavelet_map["1"] = temp_values;
186 : }
187 :
188 : // now do the cascade
189 94 : binary_vec.emplace_back("1");
190 252 : for (unsigned i=1; i<recursion_number; ++i) {
191 : std::vector<std::string> new_binary_vec;
192 1144 : for (const auto& binary : binary_vec) {
193 2958 : for (unsigned k=0; k<2; ++k) {
194 : // prepend the new bit
195 3944 : std::string new_binary = std::to_string(k) + binary;
196 1972 : mult(h_Matvec[k], scaling_map[binary], temp_values);
197 3944 : scaling_map.insert(std::pair<std::string, std::vector<double>>(new_binary, temp_values));
198 1972 : if (use_mother_wavelet) {
199 8 : mult(g_Matvec[k], scaling_map[binary], temp_values);
200 16 : wavelet_map.insert(std::pair<std::string, std::vector<double>>(new_binary, temp_values));
201 : }
202 1972 : new_binary_vec.push_back(new_binary);
203 : }
204 : }
205 158 : binary_vec = new_binary_vec;
206 158 : }
207 :
208 188 : return use_mother_wavelet ? wavelet_map : scaling_map;
209 94 : }
210 :
211 :
212 : // Fill the Grid with the values of the unordered maps
213 47 : void WaveletGrid::fillGridFromMaps(std::unique_ptr<Grid>& grid, const BinaryMap& values_map, const BinaryMap& derivs_map) {
214 47 : unsigned bins_per_int = values_map.size();
215 : // this is somewhat complicated… not sure if the unordered_map way is the best way for c++
216 1127 : for (const auto& value_iter : values_map) {
217 : // get decimal of binary key
218 1080 : unsigned decimal = std::stoi(value_iter.first, nullptr, 2);
219 : // corresponding iterator of deriv
220 1080 : auto deriv_iter = derivs_map.find(value_iter.first);
221 : // calculate first grid element
222 1080 : unsigned first_grid_element = decimal * (bins_per_int >> value_iter.first.length());
223 19480 : for (unsigned i=0; i<value_iter.second.size(); ++i) {
224 : // derivative has to be passed as vector
225 18400 : std::vector<double> deriv {deriv_iter->second[i]};
226 18400 : grid->setValueAndDerivatives(first_grid_element + bins_per_int*i, value_iter.second[i], deriv);
227 : }
228 : }
229 47 : }
230 :
231 :
232 47 : WaveletGrid::Type WaveletGrid::stringToType(std::string& type_str) {
233 188 : std::unordered_map<std::string, Type> typemap = { {"DAUBECHIES", Type::db}, {"SYMLETS", Type::sym} };
234 47 : try { return typemap.at(type_str); }
235 0 : catch(const std::out_of_range& e) {plumed_merror("The specified wavelet type "+type_str+" is not implemented."); }
236 : }
237 :
238 :
239 47 : std::string WaveletGrid::typeToString(Type type, bool abbrev) {
240 : std::string type_str;
241 47 : switch(type) {
242 23 : case Type::db:
243 23 : type_str = abbrev ? "Db" : "DAUBECHIES";
244 : break;
245 24 : case Type::sym:
246 24 : type_str = abbrev ? "Sym" : "SYMLETS";
247 : break;
248 : }
249 47 : return type_str;
250 : }
251 :
252 :
253 : }
254 : }
|
