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 "GridLinearInterpolation.h"
24 :
25 : #include "tools/Grid.h"
26 : #include "tools/Exception.h"
27 : #include "tools/Tools.h"
28 :
29 :
30 : namespace PLMD {
31 : namespace ves {
32 :
33 :
34 1869 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_1D(GridBase* grid_pntr, const std::vector<double>& arg) {
35 :
36 1869 : plumed_massert(grid_pntr->getDimension()==1,"The grid is of the wrong dimension, should be one-dimensional");
37 1869 : plumed_massert(arg.size()==1,"input value is of the wrong size");
38 :
39 1869 : double x = arg[0];
40 1869 : double grid_dx = grid_pntr->getDx()[0];
41 : double grid_min;
42 1869 : Tools::convert( grid_pntr->getMin()[0], grid_min);
43 1869 : std::vector<unsigned int> i0(1);
44 1869 : i0[0] = unsigned( std::floor( (x-grid_min)/grid_dx ) );
45 1869 : std::vector<unsigned int> i1(1);
46 1869 : i1[0] = unsigned( std::ceil( (x-grid_min)/grid_dx ) );
47 : //
48 1869 : double x0 = grid_pntr->getPoint(i0)[0];
49 1869 : double x1 = grid_pntr->getPoint(i1)[0];
50 1869 : double y0 = grid_pntr->getValue(i0);
51 1869 : double y1 = grid_pntr->getValue(i1);
52 : //
53 1869 : return linearInterpolation(x,x0,x1,y0,y1);
54 : }
55 :
56 :
57 47164 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_2D(GridBase* grid_pntr, const std::vector<double>& arg) {
58 47164 : plumed_massert(grid_pntr->getDimension()==2,"The grid is of the wrong dimension, should be two-dimensional");
59 47164 : plumed_massert(arg.size()==2,"input value is of the wrong size");
60 :
61 47164 : std::vector<double> grid_delta = grid_pntr->getDx();
62 47164 : std::vector<double> grid_min(2);
63 47164 : Tools::convert( grid_pntr->getMin()[0], grid_min[0]);
64 47164 : Tools::convert( grid_pntr->getMin()[1], grid_min[1]);
65 :
66 47164 : std::vector<unsigned int> i00(2);
67 47164 : std::vector<unsigned int> i01(2);
68 47164 : std::vector<unsigned int> i10(2);
69 47164 : std::vector<unsigned int> i11(2);
70 :
71 47164 : i00[0] = i01[0] = unsigned( std::floor( (arg[0]-grid_min[0])/grid_delta[0] ) );
72 47164 : i10[0] = i11[0] = unsigned( std::ceil( (arg[0]-grid_min[0])/grid_delta[0] ) );
73 :
74 47164 : i00[1] = i10[1] = unsigned( std::floor( (arg[1]-grid_min[1])/grid_delta[1] ) );
75 47164 : i01[1] = i11[1] = unsigned( std::ceil( (arg[1]-grid_min[1])/grid_delta[1] ) );
76 :
77 : // https://en.wikipedia.org/wiki/Bilinear_interpolation
78 47164 : double x = arg[0];
79 47164 : double y = arg[1];
80 :
81 47164 : double x0 = grid_pntr->getPoint(i00)[0];
82 47164 : double x1 = grid_pntr->getPoint(i10)[0];
83 47164 : double y0 = grid_pntr->getPoint(i00)[1];
84 47164 : double y1 = grid_pntr->getPoint(i11)[1];
85 :
86 47164 : double f00 = grid_pntr->getValue(i00);
87 47164 : double f01 = grid_pntr->getValue(i01);
88 47164 : double f10 = grid_pntr->getValue(i10);
89 47164 : double f11 = grid_pntr->getValue(i11);
90 :
91 : // linear interpolation in x-direction
92 : double fx0 = linearInterpolation(x,x0,x1,f00,f10);
93 : double fx1 = linearInterpolation(x,x0,x1,f01,f11);
94 : // linear interpolation in y-direction
95 : double fxy = linearInterpolation(y,y0,y1,fx0,fx1);
96 : //
97 47164 : return fxy;
98 : }
99 :
100 :
101 0 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_ND(GridBase* grid_pntr, const std::vector<double>& arg) {
102 0 : unsigned int dimension = grid_pntr->getDimension();
103 0 : plumed_massert(dimension==arg.size(),"The grid dimensions do not match the the given arguments.");
104 : // get points first
105 0 : std::vector<std::vector<unsigned>> point_indices = GridLinearInterpolation::getAdjacentPoints(grid_pntr, arg);
106 :
107 0 : unsigned npoints = point_indices.size();
108 : // reserve space for the point vectors and values
109 0 : std::vector<std::vector<double>> points(npoints, std::vector<double>(dimension));
110 0 : std::vector<double> values(npoints);
111 :
112 : // fill point and value vectors for all the grid points
113 0 : for (unsigned i = 0; i < npoints; ++i) {
114 0 : points[i] = grid_pntr->getPoint(point_indices[i]);
115 0 : values[i] = grid_pntr->getValue(point_indices[i]);
116 : }
117 :
118 0 : return multiLinearInterpolation(arg, points, values, dimension);
119 0 : }
120 :
121 :
122 1263854 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation_1D(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
123 1263854 : plumed_massert(grid_pntr->getDimension()==1,"The grid is of the wrong dimension, should be one-dimensional");
124 1263854 : plumed_massert(arg.size()==1,"input value is of the wrong size");
125 :
126 1263854 : double x = arg[0];
127 1263854 : double grid_dx = grid_pntr->getDx()[0];
128 : double grid_min;
129 1263854 : Tools::convert( grid_pntr->getMin()[0], grid_min);
130 :
131 1263854 : double xtoindex = (x-grid_min)/grid_dx;
132 1263854 : std::vector<unsigned int> i0(1);
133 1263854 : i0[0] = unsigned(std::floor(xtoindex));
134 1263854 : std::vector<unsigned int> i1(1);
135 1263854 : i1[0] = unsigned(std::ceil(xtoindex));
136 : //
137 1263854 : std::vector<double> d0 (1), d1 (1);
138 1263854 : double x0 = grid_pntr->getPoint(i0)[0];
139 1263854 : double x1 = grid_pntr->getPoint(i1)[0];
140 1263854 : double y0 = grid_pntr->getValueAndDerivatives(i0, d0);
141 1263854 : double y1 = grid_pntr->getValueAndDerivatives(i1, d1);
142 : //
143 1263854 : der.resize(1);
144 2500429 : der[0] = linearInterpolation(arg[0],x0,x1,d0[0],d1[0]);
145 2527708 : return linearInterpolation(arg[0],x0,x1,y0,y1);
146 : }
147 :
148 :
149 0 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation_ND(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
150 0 : unsigned int dimension = grid_pntr->getDimension();
151 0 : plumed_massert(dimension==arg.size(),"The grid dimensions do not match the given arguments");
152 : // get points first
153 0 : std::vector<std::vector<unsigned>> point_indices = getAdjacentPoints(grid_pntr, arg);
154 :
155 0 : unsigned npoints = point_indices.size();
156 : // allocate space for the point vectors and values
157 0 : std::vector<std::vector<double>> points(npoints, std::vector<double>(dimension));
158 0 : std::vector<double> values(npoints);
159 0 : std::vector<std::vector<double>> derivs(npoints, std::vector<double>(dimension));
160 :
161 : // fill point, value and deriv vectors for all the grid points
162 0 : for (unsigned i = 0; i < npoints; ++i) {
163 0 : points[i] = grid_pntr->getPoint(point_indices[i]);
164 0 : values[i] = grid_pntr->getValueAndDerivatives(point_indices[i], derivs[i]);
165 : }
166 :
167 0 : for (size_t i = 0; i < derivs.size(); ++i) {
168 0 : der[i] = multiLinearInterpolation(arg, points, derivs[i], dimension);
169 : }
170 0 : return multiLinearInterpolation(arg, points, values, dimension);
171 0 : }
172 :
173 :
174 49033 : double GridLinearInterpolation::getGridValueWithLinearInterpolation(GridBase* grid_pntr, const std::vector<double>& arg) {
175 49033 : unsigned int dim = grid_pntr->getDimension();
176 49033 : if(dim==1) {
177 1869 : return getGridValueWithLinearInterpolation_1D(grid_pntr,arg);
178 47164 : } else if(dim==2) {
179 47164 : return getGridValueWithLinearInterpolation_2D(grid_pntr,arg);
180 : } else {
181 0 : return getGridValueWithLinearInterpolation_ND(grid_pntr,arg);
182 : }
183 : }
184 :
185 :
186 1263854 : double GridLinearInterpolation::getGridValueAndDerivativesWithLinearInterpolation(GridBase* grid_pntr, const std::vector<double>& arg, std::vector<double>& der) {
187 1263854 : unsigned int dim = grid_pntr->getDimension();
188 1263854 : if(dim==1) {
189 1263854 : return getGridValueAndDerivativesWithLinearInterpolation_1D(grid_pntr,arg,der);
190 : }
191 0 : return getGridValueAndDerivativesWithLinearInterpolation_ND(grid_pntr,arg,der);
192 : }
193 :
194 :
195 : // returns the adjacent Grid indices of all double arg as vector of vectors
196 0 : std::vector<std::vector<unsigned>> GridLinearInterpolation::getAdjacentIndices(GridBase* grid_pntr, const std::vector<double>& arg) {
197 0 : unsigned int dimension = grid_pntr->getDimension();
198 :
199 0 : std::vector<std::vector<unsigned>> indices(dimension, std::vector<unsigned>(2));
200 0 : for (unsigned i=0; i<dimension; ++i) {
201 0 : std::vector<unsigned> temp_indices(2);
202 : //
203 0 : double grid_dx = grid_pntr->getDx()[i];
204 : double grid_min;
205 0 : Tools::convert( grid_pntr->getMin()[i], grid_min);
206 0 : double xtoindex = (arg[i]-grid_min)/grid_dx;
207 0 : temp_indices[0] = static_cast<unsigned>(std::floor(xtoindex));
208 0 : temp_indices[1] = static_cast<unsigned>(std::ceil(xtoindex));
209 0 : indices[i] = temp_indices;
210 : }
211 0 : return indices;
212 0 : }
213 :
214 :
215 0 : std::vector<std::vector<unsigned>> GridLinearInterpolation::getAdjacentPoints(GridBase* grid_pntr, const std::vector<double>& arg) {
216 : // upper and lower grid indices as vectors for each dimension
217 0 : std::vector<std::vector<unsigned>> grid_indices = getAdjacentIndices(grid_pntr, arg);
218 0 : unsigned npoints = 1U << grid_pntr->getDimension();
219 :
220 : // generate combination of grid indices if multidimensional to match the actual points
221 : // the retrieved combinations will be in column-major order
222 0 : std::vector<std::vector<unsigned>> point_indices(npoints);
223 0 : for (unsigned i = 0; i < npoints; ++i) {
224 : unsigned temp = i;
225 : std::vector<unsigned> current_indices;
226 0 : for (const auto& vec: grid_indices) {
227 0 : unsigned j = temp % 2;
228 0 : current_indices.push_back(vec[j]);
229 0 : temp /= 2;
230 : }
231 0 : point_indices[i] = current_indices;
232 : }
233 0 : return point_indices;
234 0 : }
235 :
236 :
237 :
238 :
239 : }
240 : }
|
