Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed 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 : plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ParallelTaskManager.h"
25 : #include "EvaluateGridFunction.h"
26 : #include "ActionWithGrid.h"
27 :
28 : //+PLUMEDOC GRIDANALYSIS INTERPOLATE_GRID
29 : /*
30 : Interpolate a smooth function stored on a grid onto a grid with a smaller grid spacing.
31 :
32 : This action takes a function evaluated on a grid as input and can be used to interpolate the values of that
33 : function on to a finer grained grid. By default the interpolation within this algorithm is done using splines.
34 :
35 : ## Examples
36 :
37 : The input below can be used to post process a trajectory. It calculates a [HISTOGRAM](HISTOGRAM.md) as a function the
38 : distance between atoms 1 and 2 using kernel density estimation. During the calculation the values of the kernels
39 : are evaluated at 100 points on a uniform grid between 0.0 and 3.0. Prior to outputting this function at the end of the
40 : simulation this function is interpolated onto a finer grid of 200 points between 0.0 and 3.0.
41 :
42 : ```plumed
43 : x: DISTANCE ATOMS=1,2
44 : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
45 : ii: INTERPOLATE_GRID ARG=hA1 GRID_BIN=200
46 : DUMPGRID ARG=ii FILE=histo.dat
47 : ```
48 :
49 : When specifying the parameters of the interpolating grid you can use GRID_BIN to specify the number of grid points
50 : or GRID_SPACING to specify the spacing between grid points as shown below:
51 :
52 : ```plumed
53 : x: DISTANCE ATOMS=1,2
54 : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
55 : ii: INTERPOLATE_GRID ARG=hA1 GRID_SPACING=0.05 INTERPOLATION_TYPE=floor
56 : DUMPGRID ARG=ii FILE=histo.dat
57 : ```
58 :
59 : In the above input spline interpolation is not used. Instead, if you evaluating the function at value $x$, which
60 : is $x_n < x < x_{n+1}$ where $x_n$ and $x_{n+1}$ are two points where the values for the function on the input grid
61 : are $f(x_n)$ and $f(x_{n+1})$ then $f(x)=f(x_n)$. If you want $f(x)=f(x_{n+1})$ you use `INTERPOLATION_TYPE=ceiling`.
62 : Alternatively, you can use linear interpolation as has been done in the following input.
63 :
64 : ```plumed
65 : x: DISTANCE ATOMS=1,2
66 : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
67 : ii: INTERPOLATE_GRID ARG=hA1 GRID_MIN=0.015 GRID_MAX=2.985 GRID_BIN=99 INTERPOLATION_TYPE=linear
68 : DUMPGRID ARG=ii FILE=histo.dat
69 : ```
70 :
71 : Notice that the way we have specified `GRID_MIN`, `GRID_MAX` and `GRID_BIN` in the above input ensures that function is evaluated
72 : at the $n-1$ mid points that lie between the points on the input grid.
73 :
74 : */
75 : //+ENDPLUMEDOC
76 :
77 : namespace PLMD {
78 : namespace gridtools {
79 :
80 : class InterpolateGrid : public ActionWithGrid {
81 : public:
82 : using input_type = EvaluateGridFunction;
83 : using PTM = ParallelTaskManager<InterpolateGrid>;
84 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
85 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
86 : private:
87 : bool firststep;
88 : /// The parallel task manager
89 : PTM taskmanager;
90 : std::vector<std::size_t> nbin;
91 : std::vector<std::string> gmin, gmax;
92 : std::vector<double> gspacing;
93 : GridCoordinatesObject output_grid;
94 : public:
95 : static void registerKeywords( Keywords& keys );
96 : explicit InterpolateGrid(const ActionOptions&ao);
97 : unsigned getNumberOfDerivatives() override ;
98 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
99 : std::vector<std::string> getGridCoordinateNames() const override ;
100 : void calculate() override ;
101 : void getArgumentInputData( std::vector<double>& inputdata ) const override ;
102 : // Calculate the value of the function at a grid point
103 : static void performTask( std::size_t task_index,
104 : const EvaluateGridFunction& actiondata,
105 : ParallelActionsInput& input,
106 : ParallelActionsOutput& output );
107 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
108 : static int getNumberOfValuesPerTask( std::size_t task_index,
109 : const EvaluateGridFunction& actiondata );
110 : static void getForceIndices( std::size_t task_index,
111 : std::size_t colno,
112 : std::size_t ntotal_force,
113 : const EvaluateGridFunction& actiondata,
114 : const ParallelActionsInput& input,
115 : ForceIndexHolder force_indices );
116 : };
117 :
118 : PLUMED_REGISTER_ACTION(InterpolateGrid,"INTERPOLATE_GRID")
119 :
120 43 : void InterpolateGrid::registerKeywords( Keywords& keys ) {
121 43 : ActionWithGrid::registerKeywords( keys );
122 43 : keys.add("optional","GRID_BIN","the number of bins for the grid");
123 86 : keys.addInputKeyword("compulsory","ARG","grid","the label for function on the grid that you would like to interpolate");
124 43 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
125 43 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid. By default this action uses the lower bound for the input grid");
126 43 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid. By default this action uses the upper bound for the input grid");
127 43 : EvaluateGridFunction::registerKeywords( keys );
128 43 : keys.remove("ZERO_OUTSIDE_GRID_RANGE");
129 86 : keys.setValueDescription("grid","the function evaluated onto the interpolated grid");
130 43 : PTM::registerKeywords( keys );
131 43 : }
132 :
133 23 : InterpolateGrid::InterpolateGrid(const ActionOptions&ao):
134 : Action(ao),
135 : ActionWithGrid(ao),
136 23 : firststep(true),
137 23 : taskmanager(this) {
138 23 : if( getNumberOfArguments()!=1 ) {
139 0 : error("should only be one argument to this action");
140 : }
141 23 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
142 0 : error("input to this action should be a grid");
143 : }
144 :
145 23 : parseVector("GRID_BIN",nbin);
146 46 : parseVector("GRID_SPACING",gspacing);
147 23 : unsigned dimension = getPntrToArgument(0)->getRank();
148 23 : gmin.resize( dimension );
149 23 : gmax.resize( dimension );
150 23 : parseVector("GRID_MIN",gmin);
151 46 : parseVector("GRID_MAX",gmax);
152 23 : if( nbin.size()!=dimension && gspacing.size()!=dimension ) {
153 0 : error("MIDPOINTS, GRID_BIN or GRID_SPACING must be set");
154 : }
155 23 : if( nbin.size()==dimension ) {
156 22 : log.printf(" number of bins in grid %zu", nbin[0]);
157 23 : for(unsigned i=1; i<nbin.size(); ++i) {
158 1 : log.printf(", %zu", nbin[i]);
159 : }
160 22 : log.printf("\n");
161 1 : } else if( gspacing.size()==dimension ) {
162 1 : log.printf(" spacing for bins in grid %f", gspacing[0]);
163 1 : for(unsigned i=1; i<gspacing.size(); ++i) {
164 0 : log.printf(", %f", gspacing[i]);
165 : }
166 1 : log.printf("\n");
167 : }
168 : // Create the input grid
169 23 : function::FunctionOptions foptions;
170 23 : EvaluateGridFunction::read( taskmanager.getActionInput(), this, foptions );
171 : // Need this for creation of tasks
172 46 : output_grid.setup( "flat", taskmanager.getActionInput().getPbc(), 0, 0.0 );
173 :
174 : // Now add a value
175 23 : std::vector<std::size_t> shape( getPntrToArgument(0)->getShape() );
176 23 : addValueWithDerivatives( shape );
177 :
178 23 : if( getPntrToArgument(0)->isPeriodic() ) {
179 : std::string min, max;
180 0 : getPntrToArgument(0)->getDomain( min, max );
181 0 : setPeriodic( min, max );
182 : } else {
183 23 : setNotPeriodic();
184 : }
185 : // Setup the task manager
186 23 : if( taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::spline ) {
187 2 : taskmanager.setupParallelTaskManager( 0, 0 );
188 21 : } else if( taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::linear ) {
189 5 : taskmanager.setupParallelTaskManager( 1+dimension, getPntrToArgument(0)->getNumberOfStoredValues() );
190 16 : } else if( taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::floor || taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::ceiling ) {
191 16 : taskmanager.setupParallelTaskManager( 1, getPntrToArgument(0)->getNumberOfStoredValues() );
192 : } else {
193 0 : error("interpolation type not defined");
194 : }
195 46 : }
196 :
197 1913 : unsigned InterpolateGrid::getNumberOfDerivatives() {
198 1913 : return getPntrToArgument(0)->getRank();
199 : }
200 :
201 343 : const GridCoordinatesObject& InterpolateGrid::getGridCoordinatesObject() const {
202 343 : return output_grid;
203 : }
204 :
205 3 : std::vector<std::string> InterpolateGrid::getGridCoordinateNames() const {
206 3 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
207 3 : plumed_assert( ag );
208 3 : return ag->getGridCoordinateNames();
209 : }
210 :
211 1763 : void InterpolateGrid::calculate() {
212 1763 : if( firststep ) {
213 23 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
214 23 : plumed_assert( ag );
215 23 : const GridCoordinatesObject& mygrid = ag->getGridCoordinatesObject();
216 47 : for(unsigned i=0; i<gmin.size(); ++i) {
217 24 : if( gmin[i]=="auto" ) {
218 36 : gmin[i] = mygrid.getMin()[i];
219 : }
220 24 : if( gmax[i]=="auto" ) {
221 36 : gmax[i] = mygrid.getMax()[i];
222 : }
223 : }
224 23 : output_grid.setBounds( gmin, gmax, nbin, gspacing );
225 23 : getPntrToComponent(0)->setShape( output_grid.getNbin(true) );
226 23 : if( taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::linear ) {
227 5 : taskmanager.setupParallelTaskManager( 1+getPntrToComponent(0)->getRank(), getPntrToArgument(0)->getNumberOfStoredValues() );
228 18 : } else if( taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::floor || taskmanager.getActionInput().interpolation_type==EvaluateGridFunction::ceiling ) {
229 16 : taskmanager.setupParallelTaskManager( 1, getPntrToArgument(0)->getNumberOfStoredValues() );
230 2 : } else if( taskmanager.getActionInput().interpolation_type!=EvaluateGridFunction::spline ) {
231 0 : error("interpolation type not defined");
232 : }
233 23 : firststep=false;
234 : }
235 1763 : taskmanager.runAllTasks();
236 1763 : }
237 :
238 1763 : void InterpolateGrid::getArgumentInputData( std::vector<double>& inputdata ) const {
239 1763 : std::size_t ndim = output_grid.getDimension();
240 1763 : std::size_t nstored = getConstPntrToComponent(0)->getNumberOfStoredValues();
241 1763 : std::vector<double> pos( ndim );
242 1763 : if( inputdata.size()!=nstored*ndim ) {
243 23 : inputdata.resize( ndim*nstored );
244 : }
245 :
246 15684 : for(unsigned i=0; i<nstored; ++i) {
247 13921 : output_grid.getGridPointCoordinates( i, pos );
248 34242 : for(unsigned j=0; j<ndim; ++j) {
249 20321 : inputdata[ i*ndim + j ] = pos[j];
250 : }
251 : }
252 1763 : }
253 :
254 15917 : void InterpolateGrid::performTask( std::size_t task_index,
255 : const EvaluateGridFunction& actiondata,
256 : ParallelActionsInput& input,
257 : ParallelActionsOutput& output ) {
258 15917 : auto funcout = function::FunctionOutput::create( 1,
259 : output.values.data(),
260 15917 : input.ranks[0],
261 : output.values.data()+1 );
262 31834 : EvaluateGridFunction::calc( actiondata,
263 : false,
264 15917 : View<const double>(input.inputdata + task_index*input.ranks[0],input.ranks[0]),
265 : funcout );
266 :
267 15917 : if( actiondata.interpolation_type==EvaluateGridFunction::linear ) {
268 6840 : View<const double> args(input.inputdata + task_index*input.ranks[0],input.ranks[0]);
269 : const GridCoordinatesObject & gridobject = actiondata.getGridObject();
270 : unsigned dimension = gridobject.getDimension();
271 6840 : std::vector<double> xfloor(dimension);
272 6840 : std::vector<unsigned> indices(dimension), nindices(dimension), ind(dimension);
273 6840 : gridobject.getIndices( args, indices );
274 6840 : unsigned nn=gridobject.getIndex(indices);
275 6840 : gridobject.getGridPointCoordinates( nn, nindices, xfloor );
276 6840 : output.derivatives[0] = 0;
277 20080 : for(unsigned i=0; i<dimension; ++i) {
278 : int x0=1;
279 13240 : if(nindices[i]==indices[i]) {
280 : x0=0;
281 : }
282 13240 : double ddx=gridobject.getGridSpacing()[i];
283 13240 : double X = fabs((args[i]-xfloor[i])/ddx-(double)x0);
284 13240 : output.derivatives[0] += (1-X);
285 13240 : output.derivatives[1+i] = X;
286 : }
287 9077 : } else if( actiondata.interpolation_type==EvaluateGridFunction::floor ) {
288 61 : output.derivatives[0] = 1;
289 9016 : } else if( actiondata.interpolation_type==EvaluateGridFunction::ceiling ) {
290 8816 : output.derivatives[0] = 1;
291 : }
292 15917 : }
293 :
294 471 : void InterpolateGrid::applyNonZeroRankForces( std::vector<double>& outforces ) {
295 471 : taskmanager.applyForces( outforces );
296 471 : }
297 :
298 1996 : int InterpolateGrid::getNumberOfValuesPerTask( std::size_t task_index,
299 : const EvaluateGridFunction& actiondata ) {
300 1996 : return 1;
301 : }
302 :
303 1996 : void InterpolateGrid::getForceIndices( std::size_t task_index,
304 : std::size_t colno,
305 : std::size_t ntotal_force,
306 : const EvaluateGridFunction& actiondata,
307 : const ParallelActionsInput& input,
308 : ForceIndexHolder force_indices ) {
309 :
310 : const GridCoordinatesObject & gridobject = actiondata.getGridObject();
311 : unsigned dimension = gridobject.getDimension();
312 1996 : std::vector<unsigned> indices(dimension);
313 1996 : gridobject.getIndices( View<const double>(input.inputdata + task_index*input.ranks[0],input.ranks[0]), indices );
314 1996 : force_indices.threadsafe_derivatives_end[0] = 0;
315 1996 : if( actiondata.interpolation_type==EvaluateGridFunction::linear ) {
316 140 : force_indices.tot_indices[0] = 1 + dimension;
317 140 : force_indices.indices[0][0] = gridobject.getIndex(indices);
318 140 : std::vector<unsigned> ind( dimension );
319 280 : for(unsigned i=0; i<dimension; ++i) {
320 280 : for(unsigned j=0; j<dimension; ++j) {
321 140 : ind[j] = indices[j];
322 : }
323 140 : if( gridobject.isPeriodic(i) && (ind[i]+1)==gridobject.getNbin(false)[i] ) {
324 0 : ind[i]=0;
325 : } else {
326 140 : ind[i] = ind[i] + 1;
327 : }
328 140 : force_indices.indices[0][1+i] = gridobject.getIndex(ind);
329 : }
330 1856 : } else if( actiondata.interpolation_type==EvaluateGridFunction::floor ) {
331 0 : force_indices.indices[0][0] = gridobject.getIndex(indices);
332 0 : force_indices.tot_indices[0] = 1;
333 1856 : } else if( actiondata.interpolation_type==EvaluateGridFunction::ceiling ) {
334 3712 : for(unsigned i=0; i<dimension; ++i) {
335 2320 : if( gridobject.isPeriodic(i) && (indices[i]+1)==gridobject.getNbin(false)[i] ) {
336 464 : indices[i]=0;
337 : } else {
338 1392 : indices[i] = indices[i] + 1;
339 : }
340 : }
341 1856 : force_indices.indices[0][0] = gridobject.getIndex(indices);
342 1856 : force_indices.tot_indices[0] = 1;
343 : } else {
344 0 : plumed_merror("cannot calculate derivatives for this type of interpolation");
345 : }
346 1996 : }
347 :
348 : }
349 : }
|
