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 "ActionWithGrid.h"
26 :
27 : //+PLUMEDOC GRIDANALYSIS MARGINAL
28 : /*
29 : Compute a marginal from the input multi-dimensional function
30 :
31 : The marginals of the function $f(x,y)$ are defined as:
32 :
33 : $$
34 : g(y) = \int f(x,y) \textrm{d}x \qquad \textrm{and} \qquad h(x) = \int f(x,y) \textrm{d}y
35 : $$
36 :
37 : This action takes a function of two or more coordinates that is stored on a grid in input. The marginal
38 : of the input function is then computed via numerical integration. To see how this works in practice
39 : consider the following input:
40 :
41 : ```plumed
42 : d1: DISTANCE ATOMS1=1,3 ATOMS2=1,4 ATOMS3=1,5 ATOMS4=1,6 COMPONENTS
43 :
44 : k: KDE ARG=d1.x,d1.y BANDWIDTH=0.02,0.02 KERNEL=gaussian GRID_BIN=50,50 GRID_MIN=0,0 GRID_MAX=1,1
45 : DUMPGRID ARG=k FILE=func_xy.grid STRIDE=1
46 : m1: MARGINAL ARG=k DIR=1
47 : DUMPGRID ARG=m1 FILE=func_x.grid STRIDE=1
48 : m2: MARGINAL ARG=k DIR=2
49 : DUMPGRID ARG=m1 FILE=func_y.grid STRIDE=1
50 : ```
51 :
52 : This input uses the [KDE](KDE.md) action to compute a distribution from the x and y components of the distances
53 : that were computed using the [DISTANCE](DISTANCE.md) command. We then output this two dimensional histogram and
54 : use the MARGINAL command to compute the distributions for the x and y components of the vector separately.
55 :
56 : */
57 : //+ENDPLUMEDOC
58 :
59 : namespace PLMD {
60 : namespace gridtools {
61 :
62 : class MarginalInput {
63 : public:
64 : /// This is the direction we are calculating the marginal along
65 : int dir=0;
66 : /// This is the function whose marginal we are calculating
67 : Value* function=NULL;
68 : /// This is the action that calculates the function we need
69 : ActionWithGrid* gridact=NULL;
70 : static void registerKeywords( Keywords& keys );
71 : static void read( MarginalInput& func, ActionWithArguments* action );
72 : // the assignment is customized (see below): copy construction is the default one
73 2 : MarginalInput() = default;
74 : MarginalInput(const MarginalInput&) = default;
75 : MarginalInput& operator=( const MarginalInput& m ) {
76 : dir = m.dir;
77 : function = m.function;
78 : gridact = m.gridact;
79 : return *this;
80 : }
81 : /// Get the vector containing the minimum value of the grid in each dimension
82 : std::vector<std::string> getMin() const ;
83 : /// Get the vector containing the maximum value of the grid in each dimension
84 : std::vector<std::string> getMax() const ;
85 : /// Get the periodicity of the grid
86 : std::vector<bool> getPbc() const ;
87 : /// Get the number of grid points in each direction
88 : std::vector<std::size_t> getNbin() const ;
89 : /// This gets the grid object
90 : const GridCoordinatesObject & getGridObject() const ;
91 : };
92 :
93 7 : void MarginalInput::registerKeywords( Keywords& keys ) {
94 7 : keys.add("compulsory","DIR","the direction along which the marginal is being computed");
95 7 : }
96 :
97 2 : std::vector<std::string> MarginalInput::getMin() const {
98 2 : std::vector<std::string> min(1, (gridact->getGridCoordinatesObject()).getMin()[dir] );
99 2 : return min;
100 : }
101 :
102 2 : std::vector<std::string> MarginalInput::getMax() const {
103 2 : std::vector<std::string> max(1, (gridact->getGridCoordinatesObject()).getMax()[dir] );
104 2 : return max;
105 : }
106 :
107 2 : std::vector<bool> MarginalInput::getPbc() const {
108 2 : std::vector<bool> pbc(1, (gridact->getGridCoordinatesObject()).isPeriodic(dir) );
109 2 : return pbc;
110 : }
111 :
112 2 : std::vector<std::size_t> MarginalInput::getNbin() const {
113 2 : std::vector<std::size_t> nbin(1, (gridact->getGridCoordinatesObject()).getNbin(false)[dir] );
114 2 : return nbin;
115 : }
116 :
117 : /// This gets the grid object
118 324 : const GridCoordinatesObject & MarginalInput::getGridObject() const {
119 324 : return (gridact->getGridCoordinatesObject());
120 : }
121 :
122 2 : void MarginalInput::read( MarginalInput& func, ActionWithArguments* action ) {
123 2 : func.function = action->getPntrToArgument(0);
124 2 : if( func.function->getRank()==0 || !func.function->hasDerivatives() ) {
125 0 : action->error("should have one grid as input to this action");
126 : }
127 2 : if( func.function->getRank()==1 ) {
128 0 : action->error("cannot calculate marginal for one dimensional grid");
129 : }
130 : // Get the input grid
131 2 : func.gridact = ActionWithGrid::getInputActionWithGrid( func.function->getPntrToAction() );
132 2 : plumed_assert( func.gridact );
133 2 : if( func.gridact->getGridCoordinatesObject().getGridType()!="flat" ) {
134 0 : action->error("cannot calculate marginals on fibonacci sphere");
135 : }
136 2 : action->parse("DIR",func.dir);
137 2 : action->log.printf(" calculating marginal along %dth axis of input grid\n", func.dir );
138 2 : func.dir = func.dir - 1;
139 2 : if( func.dir<0 || static_cast<unsigned>(func.dir)>=(action->getPntrToArgument(0))->getRank() ) {
140 0 : action->error("not enough directions in input grid to compute desired marginal");
141 : }
142 2 : }
143 :
144 : class Marginal : public ActionWithGrid {
145 : public:
146 : using input_type = MarginalInput;
147 : using PTM = ParallelTaskManager<Marginal>;
148 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
149 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
150 : private:
151 : bool firststep;
152 : /// The parallel task manager
153 : PTM taskmanager;
154 : GridCoordinatesObject output_grid;
155 : public:
156 : static void registerKeywords( Keywords& keys );
157 : explicit Marginal(const ActionOptions&ao);
158 : unsigned getNumberOfDerivatives() override ;
159 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
160 : std::vector<std::string> getGridCoordinateNames() const override ;
161 : void calculate() override ;
162 : // Calculate the value of the function at a grid point
163 : static void performTask( std::size_t task_index,
164 : const MarginalInput& actiondata,
165 : ParallelActionsInput& input,
166 : ParallelActionsOutput& output );
167 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
168 : static int getNumberOfValuesPerTask( std::size_t task_index,
169 : const MarginalInput& actiondata );
170 : static void getForceIndices( std::size_t task_index,
171 : std::size_t colno,
172 : std::size_t ntotal_force,
173 : const MarginalInput& actiondata,
174 : const ParallelActionsInput& input,
175 : ForceIndexHolder force_indices );
176 : };
177 :
178 : PLUMED_REGISTER_ACTION(Marginal,"MARGINAL")
179 :
180 7 : void Marginal::registerKeywords( Keywords& keys ) {
181 7 : ActionWithGrid::registerKeywords( keys );
182 14 : keys.addInputKeyword("compulsory","ARG","grid","the label for function on the grid that you would like to interpolate");
183 14 : keys.setValueDescription("grid","the function evaluated onto the interpolated grid");
184 7 : MarginalInput::registerKeywords( keys );
185 7 : PTM::registerKeywords( keys );
186 7 : }
187 :
188 2 : Marginal::Marginal(const ActionOptions&ao):
189 : Action(ao),
190 : ActionWithGrid(ao),
191 2 : firststep(true),
192 2 : taskmanager(this) {
193 2 : if( getNumberOfArguments()!=1 ) {
194 0 : error("should only be one argument to this action");
195 : }
196 2 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
197 0 : error("input to this action should be a grid");
198 : }
199 2 : if( getPntrToArgument(0)->getRank()==1 ) {
200 0 : error("should have more than one axis in grid to compute a marginal");
201 : }
202 :
203 : // Create the input grid
204 2 : MarginalInput::read( taskmanager.getActionInput(), this );
205 : // Need this for creation of tasks
206 4 : output_grid.setup( "flat", taskmanager.getActionInput().getPbc(), 0, 0.0 );
207 :
208 : // Now add a value
209 2 : std::vector<std::size_t> shape( 1, getPntrToArgument(0)->getShape()[taskmanager.getActionInput().dir] );
210 2 : addValueWithDerivatives( shape );
211 :
212 2 : if( getPntrToArgument(0)->isPeriodic() ) {
213 : std::string min, max;
214 0 : getPntrToArgument(0)->getDomain( min, max );
215 0 : setPeriodic( min, max );
216 : } else {
217 2 : setNotPeriodic();
218 : }
219 2 : }
220 :
221 15 : unsigned Marginal::getNumberOfDerivatives() {
222 15 : return 1;
223 : }
224 :
225 15 : const GridCoordinatesObject& Marginal::getGridCoordinatesObject() const {
226 15 : return output_grid;
227 : }
228 :
229 3 : std::vector<std::string> Marginal::getGridCoordinateNames() const {
230 3 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
231 3 : plumed_assert( ag );
232 3 : std::vector<std::string> gnames(1);
233 6 : gnames[0] = ag->getGridCoordinateNames()[taskmanager.getActionInput().dir];
234 3 : return gnames;
235 0 : }
236 :
237 6 : void Marginal::calculate() {
238 6 : if( firststep ) {
239 2 : std::vector<double> gspacing;
240 2 : output_grid.setBounds( taskmanager.getActionInput().getMin(), taskmanager.getActionInput().getMax(), taskmanager.getActionInput().getNbin(), gspacing );
241 2 : getPntrToComponent(0)->setShape( output_grid.getNbin(true) );
242 2 : taskmanager.setupParallelTaskManager( taskmanager.getActionInput().function->getNumberOfValues()/getPntrToComponent(0)->getNumberOfValues(), 0 );
243 2 : firststep=false;
244 : }
245 6 : taskmanager.runAllTasks();
246 6 : }
247 :
248 324 : void Marginal::performTask( std::size_t task_index,
249 : const MarginalInput& actiondata,
250 : ParallelActionsInput& input,
251 : ParallelActionsOutput& output ) {
252 324 : const GridCoordinatesObject & gridobject = actiondata.getGridObject();
253 324 : std::vector<unsigned> indices( gridobject.getDimension() );
254 : double volelement = 1;
255 972 : for(unsigned i=0; i<indices.size(); ++i ) {
256 648 : if( static_cast<int>(i)==actiondata.dir ) {
257 324 : continue;
258 : }
259 324 : volelement *= gridobject.getGridSpacing()[i];
260 : }
261 : std::size_t k = 0;
262 410260 : for(unsigned i=0; i<actiondata.function->getNumberOfValues(); ++i) {
263 409936 : gridobject.getIndices( i, indices );
264 409936 : if( indices[actiondata.dir]!=task_index ) {
265 403452 : continue;
266 : }
267 6484 : output.values[0] += actiondata.function->get(i);
268 6484 : output.values[1] += actiondata.function->getGridDerivative(i,0);
269 6484 : output.derivatives[k] = volelement;
270 6484 : k++;
271 : }
272 324 : output.values[0] *= volelement;
273 324 : output.values[1] *= volelement;
274 324 : }
275 :
276 0 : void Marginal::applyNonZeroRankForces( std::vector<double>& outforces ) {
277 0 : taskmanager.applyForces( outforces );
278 0 : }
279 :
280 0 : int Marginal::getNumberOfValuesPerTask( std::size_t task_index,
281 : const MarginalInput& actiondata ) {
282 0 : return 1;
283 : }
284 :
285 0 : void Marginal::getForceIndices( std::size_t task_index,
286 : std::size_t colno,
287 : std::size_t ntotal_force,
288 : const MarginalInput& actiondata,
289 : const ParallelActionsInput& input,
290 : ForceIndexHolder force_indices ) {
291 : std::size_t k = 0;
292 0 : const GridCoordinatesObject & gridobject = actiondata.getGridObject();
293 0 : std::vector<unsigned> indices( gridobject.getDimension() );
294 0 : for(unsigned i=0; i<actiondata.function->getNumberOfValues(); ++i) {
295 0 : gridobject.getIndices( i, indices );
296 0 : if( indices[actiondata.dir]!=task_index ) {
297 0 : continue;
298 : }
299 0 : force_indices.indices[0][k] = i;
300 0 : k++;
301 : }
302 0 : force_indices.threadsafe_derivatives_end[0] = k;
303 0 : force_indices.tot_indices[0] = k;
304 0 : }
305 :
306 : }
307 : }
|
