Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2020 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 : #ifndef __PLUMED_function_FunctionOfVector_h
23 : #define __PLUMED_function_FunctionOfVector_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "core/ActionSetup.h"
28 : #include "FunctionSetup.h"
29 : #include "Custom.h"
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 : template <class CV, typename myPTM=defaultPTM>
35 : class FunctionOfVector : public ActionWithVector {
36 : public:
37 : using input_type = FunctionData<CV>;
38 : using mytype= FunctionOfVector<CV,myPTM>;
39 : using PTM = typename myPTM::template PTM<mytype>;
40 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
41 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
42 : private:
43 : /// The parallel task manager
44 : PTM taskmanager;
45 : /// Set equal to one if we are doing EvaluateGridFunction
46 : unsigned argstart;
47 : /// Can an element be skipped when all its arguments are zero (i.e. is the function zero there)?
48 : bool skipZeroElements;
49 : public:
50 : static void registerKeywords(Keywords&);
51 : explicit FunctionOfVector(const ActionOptions&);
52 3258 : ~FunctionOfVector() {}
53 : std::string getOutputComponentDescription( const std::string& cname,
54 : const Keywords& keys ) const override ;
55 : /// Get the number of derivatives for this action
56 : unsigned getNumberOfDerivatives() override ;
57 : /// Resize vectors that are the wrong size
58 : void prepare() override ;
59 : /// Get the label to write in the graph
60 : std::string writeInGraph() const override ;
61 : /// This builds the task list for the action
62 : void calculate() override;
63 : /// Elements whose arguments are all zero are skipped only if the function is zero there
64 : int checkTaskIsActive( const unsigned& itask ) const override ;
65 : /// Add some forces
66 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
67 : /// Get the input data
68 : void getArgumentInputData( std::vector<double>& inputdata ) const override ;
69 : void getArgumentInputData( std::vector<float>& inputdata ) const override ;
70 : /// Calculate the function
71 : static void performTask( std::size_t task_index,
72 : const FunctionData<CV>& actiondata,
73 : ParallelActionsInput& input,
74 : ParallelActionsOutput& output );
75 : /// Get the indices of the forces
76 : static int getNumberOfValuesPerTask( std::size_t task_index, const FunctionData<CV>& actiondata );
77 : static void getForceIndices( std::size_t task_index,
78 : std::size_t colno,
79 : std::size_t ntotal_force,
80 : const FunctionData<CV>& actiondata,
81 : const ParallelActionsInput& input,
82 : ForceIndexHolder force_indices );
83 : };
84 :
85 : template <class CV, typename myPTM>
86 3317 : void FunctionOfVector<CV, myPTM>::registerKeywords(Keywords& keys ) {
87 3317 : Action::registerKeywords(keys);
88 3317 : ActionWithValue::registerKeywords(keys);
89 3317 : ActionWithArguments::registerKeywords(keys);
90 3317 : std::string name = keys.getDisplayName();
91 3317 : std::size_t und=name.find("_VECTOR");
92 3317 : keys.setDisplayName( name.substr(0,und) );
93 6634 : keys.addInputKeyword("compulsory","ARG","scalar/vector","the labels of the scalar and vector that on which the function is being calculated elementwise");
94 3317 : keys.reserve("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function. If the output is not periodic you must state this using PERIODIC=NO");
95 3317 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
96 3317 : CV::registerKeywords( keys );
97 3317 : if( keys.getDisplayName()=="SUM" ) {
98 0 : keys.setValueDescription("scalar","the sum of all the elements in the input vector");
99 3317 : } else if( keys.getDisplayName()=="MEAN" ) {
100 0 : keys.setValueDescription("scalar","the mean of all the elements in the input vector");
101 3317 : } else if( keys.getDisplayName()=="HIGHEST" ) {
102 8 : keys.setValueDescription("scalar/vector","the largest element of the input vector if one vector specified. If multiple vectors of the same size specified the largest elements of these vector computed elementwise.");
103 8 : keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
104 3313 : } else if( keys.getDisplayName()=="LOWEST" ) {
105 116 : keys.setValueDescription("scalar/vector","the smallest element in the input vector if one vector specified. If multiple vectors of the same size specified the smallest elements of these vector computed elementwise.");
106 116 : keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
107 3255 : } else if( keys.getDisplayName()=="SORT" ) {
108 14 : keys.setValueDescription("vector","a vector that has been sorted into ascending order");
109 14 : keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
110 6496 : } else if( keys.outputComponentExists(".#!value") ) {
111 6488 : keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed");
112 12976 : keys.setValueDescription("vector","the vector obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input vectors");
113 : }
114 3317 : if( keys.getDisplayName()=="EVALUATE_FUNCTION_FROM_GRID" ) {
115 16 : keys.addInputKeyword("compulsory","ARG","scalar/vector/grid","");
116 : }
117 3317 : PTM::registerKeywords( keys );
118 3317 : }
119 :
120 : template <class CV, typename myPTM>
121 1629 : FunctionOfVector<CV, myPTM>::FunctionOfVector(const ActionOptions&ao):
122 : Action(ao),
123 : ActionWithVector(ao),
124 1629 : taskmanager(this),
125 1629 : argstart(0),
126 1629 : skipZeroElements(true) {
127 : // Check if first argument is grid
128 1629 : if( getPntrToArgument(0)->getRank()>0 && getPntrToArgument(0)->hasDerivatives() ) {
129 2 : argstart=1;
130 : }
131 1629 : if( getNumberOfArguments()==argstart ) {
132 0 : error("no arguments specified");
133 : }
134 :
135 1629 : if( getPntrToArgument(argstart)->getRank()!=1 ) {
136 0 : error("first argument to this action must be a vector");
137 : }
138 :
139 : // Get the number of arguments
140 : unsigned nargs = getNumberOfArguments();
141 : int nmasks = getNumberOfMasks();
142 1629 : if( nargs>=static_cast<unsigned>(nmasks) && nmasks>0 ) {
143 91 : nargs = nargs - nmasks;
144 : }
145 :
146 : // Get the shape of the output
147 : std::size_t nscalars = 0;
148 1629 : std::vector<std::size_t> shape(1);
149 1629 : shape[0]=getPntrToArgument(argstart)->getShape()[0];
150 3498 : for(unsigned i=argstart+1; i<nargs; ++i) {
151 1869 : if( getPntrToArgument(i)->getRank()==0 ) {
152 247 : nscalars++;
153 1622 : } else if( getPntrToArgument(i)->getRank()==1 ) {
154 1622 : if( getPntrToArgument(i)->getShape()[0]!=shape[0] ) {
155 0 : error("mismatch between sizes of input arguments");
156 1622 : } else if( nscalars>0 ) {
157 0 : error("scalars should be specified in argument list after all vectors");
158 : }
159 : } else {
160 0 : error("input arguments should be vectors or scalars");
161 : }
162 : }
163 1629 : if( nmasks>0 && getPntrToArgument(getNumberOfArguments()-nmasks)->getShape()[0]!=shape[0] ) {
164 0 : error("input mask has wrong size");
165 : }
166 :
167 : // Setup the function and the values
168 : // Get the function data from the parallel task manager, to avoid copies
169 1629 : auto & myfunc = taskmanager.getActionInput();
170 1629 : myfunc.argstart = argstart;
171 1629 : myfunc.nscalars = nscalars;
172 1629 : FunctionData<CV>::setup( myfunc.f, keywords.getOutputComponents(), shape, false, this );
173 : // ActionWithVector skips the elements for which all the arguments are zero (and have zero
174 : // derivatives there); the output is then zero, which is wrong if the function is not zero
175 : // there (e.g. FUNC=x+1). The check is only needed when no masks, scalars or grids are used:
176 : // in the other cases the elements are not skipped this way.
177 1629 : if( nmasks<=0 && argstart==0 && nscalars==0 ) {
178 1290 : skipZeroElements = functionIsZeroWhenArgumentsAreZero( myfunc.f, nargs, getNumberOfComponents() );
179 1290 : if( !skipZeroElements ) {
180 212 : log.printf(" the function is not zero when its arguments are zero: all the elements are computed\n");
181 : }
182 : }
183 : // Setup the parallel task manager
184 1629 : taskmanager.setupParallelTaskManager( nargs-argstart, nscalars );
185 1629 : }
186 :
187 : template <class CV, typename myPTM>
188 11836970 : int FunctionOfVector<CV, myPTM>::checkTaskIsActive( const unsigned& itask ) const {
189 11836970 : if( !skipZeroElements ) {
190 : return 1;
191 : }
192 10402778 : return ActionWithVector::checkTaskIsActive( itask );
193 : }
194 :
195 : template <class CV, typename myPTM>
196 2 : std::string FunctionOfVector<CV,
197 : myPTM>::getOutputComponentDescription( const std::string& cname,
198 : const Keywords& keys ) const {
199 2 : if( getName().find("SORT")==std::string::npos ) {
200 0 : return ActionWithValue::getOutputComponentDescription( cname, keys );
201 : }
202 4 : return "the " + cname + "th smallest element in the input vectors";
203 : }
204 :
205 : template <class CV, typename myPTM>
206 11 : std::string FunctionOfVector<CV, myPTM>::writeInGraph() const {
207 11 : std::size_t und = getName().find_last_of("_");
208 11 : return getName().substr(0,und);
209 : }
210 :
211 : template <class CV, typename myPTM>
212 5869 : unsigned FunctionOfVector<CV, myPTM>::getNumberOfDerivatives() {
213 : unsigned nder = 0;
214 20004 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
215 14135 : nder += getPntrToArgument(i)->getNumberOfStoredValues();
216 : }
217 5869 : return nder;
218 : }
219 :
220 : template <class CV, typename myPTM>
221 103132 : void FunctionOfVector<CV, myPTM>::prepare() {
222 103132 : std::vector<std::size_t> shape(1);
223 103132 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
224 103132 : if( getPntrToArgument(i)->getRank()==1 ) {
225 103132 : shape[0] = getPntrToArgument(i)->getShape()[0];
226 103132 : break;
227 : }
228 : }
229 206264 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
230 103132 : Value* myval = getPntrToComponent(i);
231 103132 : if( myval->getRank()==1 && myval->getShape()[0]!=shape[0] ) {
232 47 : myval->setShape(shape);
233 : }
234 : }
235 103132 : ActionWithVector::prepare();
236 103132 : }
237 :
238 : template <class CV, typename myPTM>
239 102523 : void FunctionOfVector<CV, myPTM>::getArgumentInputData( std::vector<double>& inputdata ) const {
240 102523 : unsigned nargs = getNumberOfArguments();
241 : int nmasks = getNumberOfMasks();
242 102523 : if( nargs>=static_cast<unsigned>(nmasks) && nmasks>0 ) {
243 554 : nargs = nargs - nmasks;
244 : }
245 :
246 : std::size_t ntasks = 0;
247 102523 : for(unsigned i=argstart; i<nargs; ++i) {
248 102523 : if( getPntrToArgument(i)->getRank()==1 ) {
249 102523 : ntasks = getPntrToArgument(i)->getShape()[0];
250 102523 : break;
251 : }
252 : }
253 :
254 102523 : std::size_t ndata = static_cast<std::size_t>(nargs-argstart)*ntasks;
255 102523 : if( inputdata.size()!=ndata ) {
256 1499 : inputdata.resize( ndata );
257 : }
258 :
259 278141 : for(unsigned j=argstart; j<nargs; ++j) {
260 : const Value* myarg = getPntrToArgument(j);
261 175618 : if( myarg->getRank()==0 ) {
262 54372 : double val = myarg->get();
263 2314877 : for(unsigned i=0; i<ntasks; ++i) {
264 2260505 : inputdata[(nargs-argstart)*i + j-argstart] = val;
265 : }
266 : } else {
267 23076786 : for(unsigned i=0; i<ntasks; ++i) {
268 22955540 : inputdata[(nargs-argstart)*i + j-argstart] = myarg->get(i);
269 : }
270 : }
271 : }
272 102523 : }
273 :
274 : template <class CV, typename myPTM>
275 0 : void FunctionOfVector<CV, myPTM>::getArgumentInputData( std::vector<float>& inputdata ) const {
276 0 : unsigned nargs = getNumberOfArguments();
277 : int nmasks = getNumberOfMasks();
278 0 : if( nargs>=static_cast<unsigned>(nmasks) && nmasks>0 ) {
279 0 : nargs = nargs - nmasks;
280 : }
281 :
282 : std::size_t ntasks = 0;
283 0 : for(unsigned i=argstart; i<nargs; ++i) {
284 0 : if( getPntrToArgument(i)->getRank()==1 ) {
285 0 : ntasks = getPntrToArgument(i)->getShape()[0];
286 0 : break;
287 : }
288 : }
289 :
290 0 : std::size_t ndata = static_cast<std::size_t>(nargs-argstart)*ntasks;
291 0 : if( inputdata.size()!=ndata ) {
292 0 : inputdata.resize( ndata );
293 : }
294 :
295 0 : for(unsigned j=argstart; j<nargs; ++j) {
296 : const Value* myarg = getPntrToArgument(j);
297 0 : if( myarg->getRank()==0 ) {
298 0 : double val = myarg->get();
299 0 : for(unsigned i=0; i<ntasks; ++i) {
300 0 : inputdata[(nargs-argstart)*i + j-argstart] = val;
301 : }
302 : } else {
303 0 : for(unsigned i=0; i<ntasks; ++i) {
304 0 : inputdata[(nargs-argstart)*i + j-argstart] = myarg->get(i);
305 : }
306 : }
307 : }
308 0 : }
309 :
310 : template <class CV, typename myPTM>
311 102523 : void FunctionOfVector<CV, myPTM>::calculate() {
312 : // This is done if we are calculating a function of multiple cvs
313 102523 : taskmanager.runAllTasks();
314 102523 : }
315 :
316 : template <class CV, typename myPTM>
317 14323473 : void FunctionOfVector<CV, myPTM>::performTask( std::size_t task_index,
318 : const FunctionData<CV>& actiondata,
319 : ParallelActionsInput& input,
320 : ParallelActionsOutput& output ) {
321 14323473 : auto funcout = FunctionOutput::create( input.ncomponents,
322 : output.values.data(),
323 14323473 : input.nderivatives_per_scalar,
324 : output.derivatives.data() );
325 14341388 : CV::calc( actiondata.f,
326 14323473 : input.noderiv,
327 14323473 : View<const double>( input.inputdata + task_index*input.nderivatives_per_scalar,
328 : input.nderivatives_per_scalar ),
329 : funcout );
330 14323473 : }
331 :
332 : template <class CV, typename myPTM>
333 88009 : void FunctionOfVector<CV, myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
334 88009 : taskmanager.applyForces( outforces );
335 88009 : }
336 :
337 : template <class CV, typename myPTM>
338 5405210 : int FunctionOfVector<CV,
339 : myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
340 : const FunctionData<CV>& actiondata ) {
341 5405210 : return 1;
342 : }
343 :
344 : template <class CV, typename myPTM>
345 7967388 : void FunctionOfVector<CV, myPTM>::getForceIndices( std::size_t task_index,
346 : std::size_t colno,
347 : std::size_t ntotal_force,
348 : const FunctionData<CV>& actiondata,
349 : const ParallelActionsInput& input,
350 : ForceIndexHolder force_indices ) {
351 :
352 7967388 : unsigned vector_end = actiondata.argstart + input.nderivatives_per_scalar - actiondata.nscalars;
353 23314784 : for(unsigned j=0; j<input.ncomponents; ++j) {
354 60386930 : for(unsigned k=actiondata.argstart; k<vector_end; ++k) {
355 45039534 : force_indices.indices[j][k-actiondata.argstart] = input.argstarts[k] + task_index;
356 : }
357 17424884 : for(unsigned k=vector_end; k<vector_end+actiondata.nscalars; ++k) {
358 2077488 : force_indices.indices[j][k-actiondata.argstart] = input.argstarts[k];
359 : }
360 15347396 : force_indices.threadsafe_derivatives_end[j] = input.nderivatives_per_scalar-actiondata.nscalars;
361 15347396 : force_indices.tot_indices[j] = input.nderivatives_per_scalar;
362 : }
363 7967388 : }
364 :
365 :
366 :
367 : /// Explicit specialization defined in the corresponding .cpp file; declared
368 : /// here so that every translation unit instantiating the template sees it.
369 : template <> std::string FunctionOfVector<Custom>::writeInGraph() const;
370 : }
371 : }
372 : #endif
|
