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/CollectFrames.h"
27 : #include "core/ActionSetup.h"
28 : #include "tools/Matrix.h"
29 : #include "Sum.h"
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 : template <class T>
35 : class FunctionOfVector : public ActionWithVector {
36 : private:
37 : /// Do the calculation at the end of the run
38 : bool doAtEnd;
39 : /// Is this the first time we are doing the calc
40 : bool firststep;
41 : /// The function that is being computed
42 : T myfunc;
43 : /// The number of derivatives for this action
44 : unsigned nderivatives;
45 : /// A vector that tells us if we have stored the input value
46 : std::vector<bool> stored_arguments;
47 : public:
48 : static void registerKeywords(Keywords&);
49 : /// This method is used to run the calculation with functions such as highest/lowest and sort.
50 : /// It is static so we can reuse the functionality in FunctionOfMatrix
51 : static void runSingleTaskCalculation( const Value* arg, ActionWithValue* action, T& f );
52 : explicit FunctionOfVector(const ActionOptions&);
53 4752 : ~FunctionOfVector() {}
54 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
55 : /// Get the size of the task list at the end of the run
56 : unsigned getNumberOfFinalTasks();
57 : /// Check if derivatives are available
58 : void turnOnDerivatives() override;
59 : /// Get the number of derivatives for this action
60 : unsigned getNumberOfDerivatives() override ;
61 : /// Resize vectors that are the wrong size
62 : void prepare() override ;
63 : /// Check if all he actions are required
64 : void areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions );
65 : /// Get the label to write in the graph
66 20 : std::string writeInGraph() const override {
67 20 : return myfunc.getGraphInfo( getName() );
68 : }
69 : /// This builds the task list for the action
70 : void calculate() override;
71 : /// This ensures that we create some bookeeping stuff during the first step
72 : void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
73 : /// Calculate the function
74 : void performTask( const unsigned& current, MultiValue& myvals ) const override ;
75 : };
76 :
77 : template <class T>
78 4690 : void FunctionOfVector<T>::registerKeywords(Keywords& keys ) {
79 4690 : Action::registerKeywords(keys);
80 4690 : ActionWithValue::registerKeywords(keys);
81 4690 : ActionWithArguments::registerKeywords(keys);
82 4690 : keys.use("ARG");
83 4690 : std::string name = keys.getDisplayName();
84 4690 : std::size_t und=name.find("_VECTOR");
85 4690 : keys.setDisplayName( name.substr(0,und) );
86 9380 : 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");
87 9380 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
88 2081 : T tfunc;
89 4690 : tfunc.registerKeywords( keys );
90 9380 : if( keys.getDisplayName()=="SUM" ) {
91 2480 : keys.setValueDescription("the sum of all the elements in the input vector");
92 6900 : } else if( keys.getDisplayName()=="MEAN" ) {
93 718 : keys.setValueDescription("the mean of all the elements in the input vector");
94 6182 : } else if( keys.getDisplayName()=="HIGHEST" ) {
95 88 : keys.setValueDescription("the largest element of the input vector");
96 6094 : } else if( keys.getDisplayName()=="LOWEST" ) {
97 134 : keys.setValueDescription("the smallest element in the input vector");
98 5960 : } else if( keys.getDisplayName()=="SORT" ) {
99 28 : keys.setValueDescription("a vector that has been sorted into ascending order");
100 5932 : } else if( keys.outputComponentExists(".#!value") ) {
101 5908 : keys.setValueDescription("the vector obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input vectors");
102 : }
103 7240 : }
104 :
105 : template <class T>
106 2307 : FunctionOfVector<T>::FunctionOfVector(const ActionOptions&ao):
107 : Action(ao),
108 : ActionWithVector(ao),
109 2307 : doAtEnd(true),
110 2307 : firststep(true),
111 2307 : nderivatives(0) {
112 : // Get the shape of the output
113 2307 : std::vector<unsigned> shape(1);
114 2307 : shape[0]=getNumberOfFinalTasks();
115 : // Read the input and do some checks
116 2307 : myfunc.read( this );
117 : // Create the task list
118 2164 : if( myfunc.doWithTasks() ) {
119 2267 : doAtEnd=false;
120 2267 : if( shape[0]>0 ) {
121 2204 : done_in_chain=true;
122 : }
123 : } else {
124 40 : plumed_assert( getNumberOfArguments()==1 );
125 40 : done_in_chain=false;
126 40 : getPntrToArgument(0)->buildDataStore();
127 : }
128 : // Get the names of the components
129 2307 : std::vector<std::string> components( keywords.getOutputComponents() );
130 : // Create the values to hold the output
131 59 : std::vector<std::string> str_ind( myfunc.getComponentsPerLabel() );
132 4614 : for(unsigned i=0; i<components.size(); ++i) {
133 2307 : if( str_ind.size()>0 ) {
134 8 : std::string strcompn = components[i];
135 8 : if( components[i]==".#!value" ) {
136 : strcompn = "";
137 : }
138 34 : for(unsigned j=0; j<str_ind.size(); ++j) {
139 7 : if( myfunc.zeroRank() ) {
140 52 : addComponentWithDerivatives( strcompn + str_ind[j] );
141 : } else {
142 0 : addComponent( strcompn + str_ind[j], shape );
143 : }
144 : }
145 2299 : } else if( components[i].find_first_of("_")!=std::string::npos ) {
146 0 : if( getNumberOfArguments()==1 && myfunc.zeroRank() ) {
147 0 : addValueWithDerivatives();
148 0 : } else if( getNumberOfArguments()==1 ) {
149 0 : addValue( shape );
150 : } else {
151 : unsigned argstart=myfunc.getArgStart();
152 0 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
153 0 : if( myfunc.zeroRank() ) {
154 0 : addComponentWithDerivatives( getPntrToArgument(i)->getName() + components[i] );
155 : } else {
156 0 : addComponent( getPntrToArgument(i)->getName() + components[i], shape );
157 : }
158 : }
159 : }
160 845 : } else if( components[i]==".#!value" && myfunc.zeroRank() ) {
161 1654 : addValueWithDerivatives();
162 1472 : } else if( components[i]==".#!value" ) {
163 1472 : addValue(shape);
164 0 : } else if( myfunc.zeroRank() ) {
165 0 : addComponentWithDerivatives( components[i] );
166 : } else {
167 0 : addComponent( components[i], shape );
168 : }
169 : }
170 : // Check if we can turn off the derivatives when they are zero
171 1033 : if( myfunc.getDerivativeZeroIfValueIsZero() ) {
172 624 : for(int i=0; i<getNumberOfComponents(); ++i) {
173 312 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
174 : }
175 : }
176 : // Check if this is a timeseries
177 : unsigned argstart=myfunc.getArgStart();
178 : // for(unsigned i=argstart; i<getNumberOfArguments();++i) {
179 : // if( getPntrToArgument(i)->isTimeSeries() ) {
180 : // for(unsigned i=0; i<getNumberOfComponents(); ++i) getPntrToOutput(i)->makeHistoryDependent();
181 : // break;
182 : // }
183 : // }
184 : // Set the periodicities of the output components
185 2307 : myfunc.setPeriodicityForOutputs( this );
186 : // Check if we can put the function in a chain
187 6235 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
188 : // CollectFrames* ab=dynamic_cast<CollectFrames*>( getPntrToArgument(i)->getPntrToAction() );
189 : // if( ab && ab->hasClear() ) { doNotChain=true; getPntrToArgument(i)->buildDataStore( getLabel() ); }
190 : // No chains if we are using a sum or a mean
191 3928 : if( getPntrToArgument(i)->getRank()==0 ) {
192 246 : FunctionOfVector<Sum>* as = dynamic_cast<FunctionOfVector<Sum>*>( getPntrToArgument(i)->getPntrToAction() );
193 246 : if(as) {
194 48 : done_in_chain=false;
195 : }
196 : } else {
197 3682 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
198 3682 : if( !av ) {
199 450 : done_in_chain=false;
200 : }
201 : }
202 : }
203 : unsigned ignore, atstart = argstart;
204 : ActionAtomistic* firstat = NULL;
205 4383 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
206 2847 : ActionWithVector* myvec = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
207 2847 : if( !myvec ) {
208 623 : continue;
209 : }
210 2224 : firstat = dynamic_cast<ActionAtomistic*>( myvec->getFirstActionInChain() );
211 2224 : if( firstat && myvec->isInSubChain(ignore) ) {
212 : break;
213 : }
214 : atstart = i;
215 : }
216 5964 : for(unsigned i=atstart; i<getNumberOfArguments(); ++i) {
217 3697 : ActionWithVector* myvec = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
218 3697 : if( !myvec ) {
219 596 : continue;
220 : }
221 3101 : ActionAtomistic* secondat = dynamic_cast<ActionAtomistic*>( myvec->getFirstActionInChain() );
222 3101 : if( secondat!=firstat && myvec->isInSubChain(ignore) ) {
223 40 : done_in_chain = false;
224 40 : break;
225 : }
226 : }
227 : // Don't need to do the calculation in a chain if the input is constant
228 : bool allconstant=true;
229 2646 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
230 2318 : if( !getPntrToArgument(i)->isConstant() ) {
231 : allconstant=false;
232 : break;
233 : }
234 : }
235 2307 : if( allconstant ) {
236 328 : done_in_chain=false;
237 : }
238 2307 : nderivatives = buildArgumentStore(myfunc.getArgStart());
239 4614 : }
240 :
241 : template <class T>
242 5 : std::string FunctionOfVector<T>::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
243 5 : if( getName().find("SORT")==std::string::npos ) {
244 0 : return ActionWithValue::getOutputComponentDescription( cname, keys );
245 : }
246 5 : if( getNumberOfArguments()==1 ) {
247 6 : return "the " + cname + "th largest element of the vector " + getPntrToArgument(0)->getName();
248 : }
249 4 : return "the " + cname + "th largest element in the input vectors";
250 : }
251 :
252 : template <class T>
253 6208 : void FunctionOfVector<T>::turnOnDerivatives() {
254 6208 : if( !getPntrToComponent(0)->isConstant() && !myfunc.derivativesImplemented() ) {
255 0 : error("derivatives have not been implemended for " + getName() );
256 : }
257 6208 : ActionWithValue::turnOnDerivatives();
258 6208 : myfunc.setup(this );
259 6208 : }
260 :
261 : template <class T>
262 53246 : unsigned FunctionOfVector<T>::getNumberOfDerivatives() {
263 53246 : return nderivatives;
264 : }
265 :
266 : template <class T>
267 186141 : void FunctionOfVector<T>::prepare() {
268 : unsigned argstart = myfunc.getArgStart();
269 186141 : std::vector<unsigned> shape(1);
270 233119 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
271 233119 : if( getPntrToArgument(i)->getRank()==1 ) {
272 186141 : shape[0] = getPntrToArgument(i)->getShape()[0];
273 186141 : break;
274 : }
275 : }
276 372698 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
277 186557 : Value* myval = getPntrToComponent(i);
278 186557 : if( myval->getRank()==1 && myval->getShape()[0]!=shape[0] ) {
279 47 : myval->setShape(shape);
280 : }
281 : }
282 186141 : ActionWithVector::prepare();
283 186141 : }
284 :
285 : template <class T>
286 313237 : void FunctionOfVector<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
287 313237 : if( firststep ) {
288 2199 : stored_arguments.resize( getNumberOfArguments() );
289 2199 : std::string control = getFirstActionInChain()->getLabel();
290 6000 : for(unsigned i=0; i<stored_arguments.size(); ++i) {
291 3801 : if( getPntrToArgument(i)->isConstant() ) {
292 : stored_arguments[i]=false;
293 : } else {
294 3274 : stored_arguments[i] = !getPntrToArgument(i)->ignoreStoredValue( control );
295 : }
296 : }
297 2199 : firststep=false;
298 : }
299 313237 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
300 313237 : }
301 :
302 : template <class T>
303 7190729 : void FunctionOfVector<T>::performTask( const unsigned& current, MultiValue& myvals ) const {
304 8375 : unsigned argstart=myfunc.getArgStart();
305 7190729 : std::vector<double> args( getNumberOfArguments()-argstart);
306 7190729 : if( actionInChain() ) {
307 9638352 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
308 6362122 : if( getPntrToArgument(i)->getRank()==0 ) {
309 17788 : args[i-argstart] = getPntrToArgument(i)->get();
310 6344334 : } else if( !getPntrToArgument(i)->valueHasBeenSet() ) {
311 6187993 : args[i-argstart] = myvals.get( getPntrToArgument(i)->getPositionInStream() );
312 : } else {
313 156341 : args[i-argstart] = getPntrToArgument(i)->get( myvals.getTaskIndex() );
314 : }
315 : }
316 : } else {
317 11599144 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
318 7684645 : if( getPntrToArgument(i)->getRank()==1 ) {
319 4442100 : args[i-argstart]=getPntrToArgument(i)->get(current);
320 : } else {
321 3242545 : args[i-argstart] = getPntrToArgument(i)->get();
322 : }
323 : }
324 : }
325 : // Calculate the function and its derivatives
326 7190729 : std::vector<double> vals( getNumberOfComponents() );
327 7190729 : Matrix<double> derivatives( getNumberOfComponents(), args.size() );
328 7190729 : myfunc.calc( this, args, vals, derivatives );
329 : // And set the values
330 14381458 : for(unsigned i=0; i<vals.size(); ++i) {
331 7190729 : myvals.addValue( getConstPntrToComponent(i)->getPositionInStream(), vals[i] );
332 : }
333 : // Return if we are not computing derivatives
334 7190729 : if( doNotCalculateDerivatives() ) {
335 : return;
336 : }
337 : // And now compute the derivatives
338 : // Second condition here is only not true if actionInChain()==True if
339 : // input arguments the only non-constant objects in input are scalars.
340 : // In that case we can use the non chain version to calculate the derivatives
341 : // with respect to the scalar.
342 6274299 : if( actionInChain() ) {
343 7990707 : for(unsigned j=0; j<args.size(); ++j) {
344 8375 : unsigned istrn = getPntrToArgument(argstart+j)->getPositionInStream();
345 5326473 : if( stored_arguments[argstart+j] ) {
346 70 : unsigned task_index = myvals.getTaskIndex();
347 70 : if( getPntrToArgument(argstart+j)->getRank()==0 ) {
348 : task_index=0;
349 : }
350 70 : myvals.addDerivative( istrn, task_index, 1.0 );
351 70 : myvals.updateIndex( istrn, task_index );
352 : }
353 5326473 : unsigned arg_deriv_s = arg_deriv_starts[argstart+j];
354 173086446 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
355 167759973 : unsigned kind=myvals.getActiveIndex(istrn,k);
356 335519946 : for(int i=0; i<getNumberOfComponents(); ++i) {
357 167759973 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
358 167759973 : myvals.addDerivative( ostrn, arg_deriv_s + kind, derivatives(i,j)*myvals.getDerivative( istrn, kind ) );
359 : }
360 : }
361 : // Ensure we only store one lot of derivative indices
362 : bool found=false;
363 5326473 : ActionWithValue* aav=getPntrToArgument(argstart+j)->getPntrToAction();
364 5347546 : for(unsigned k=0; k<j; ++k) {
365 2662303 : if( arg_deriv_starts[argstart+k]==arg_deriv_s ) {
366 2641230 : if( getPntrToArgument(argstart+k)->getPntrToAction()!=aav ) {
367 386454 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(argstart+j)->getPntrToAction() );
368 386454 : if( av ) {
369 772908 : for(int i=0; i<getNumberOfComponents(); ++i) {
370 386454 : av->updateAdditionalIndices( getConstPntrToComponent(i)->getPositionInStream(), myvals );
371 : }
372 : }
373 : }
374 : found=true;
375 : break;
376 : }
377 : }
378 2641230 : if( found ) {
379 2641230 : continue;
380 : }
381 88267362 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
382 85582119 : unsigned kind=myvals.getActiveIndex(istrn,k);
383 171164238 : for(int i=0; i<getNumberOfComponents(); ++i) {
384 85582119 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
385 85582119 : myvals.updateIndex( ostrn, arg_deriv_s + kind );
386 : }
387 : }
388 : }
389 : } else {
390 : unsigned base=0;
391 10697896 : for(unsigned j=0; j<args.size(); ++j) {
392 7087831 : if( getPntrToArgument(argstart+j)->getRank()==1 ) {
393 7983052 : for(int i=0; i<getNumberOfComponents(); ++i) {
394 3991526 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
395 3991526 : myvals.addDerivative( ostrn, base+current, derivatives(i,j) );
396 3991526 : myvals.updateIndex( ostrn, base+current );
397 : }
398 : } else {
399 6192610 : for(int i=0; i<getNumberOfComponents(); ++i) {
400 3096305 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
401 3096305 : myvals.addDerivative( ostrn, base, derivatives(i,j) );
402 3096305 : myvals.updateIndex( ostrn, base );
403 : }
404 : }
405 7087831 : base += getPntrToArgument(argstart+j)->getNumberOfValues();
406 : }
407 : }
408 : }
409 :
410 : template <class T>
411 2307 : unsigned FunctionOfVector<T>::getNumberOfFinalTasks() {
412 : unsigned nelements=0, argstart=myfunc.getArgStart();
413 6235 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
414 3928 : plumed_assert( getPntrToArgument(i)->getRank()<2 );
415 3928 : if( getPntrToArgument(i)->getRank()==1 ) {
416 3682 : if( nelements>0 ) {
417 : // if( getPntrToArgument(i)->isTimeSeries() && getPntrToArgument(i)->getShape()[0]<nelements ) nelements=getPntrToArgument(i)->isTimeSeries();
418 : // else
419 1375 : if(getPntrToArgument(i)->getShape()[0]!=nelements ) {
420 0 : error("all vectors input should have the same length");
421 : }
422 : } else if( nelements==0 ) {
423 2307 : nelements=getPntrToArgument(i)->getShape()[0];
424 : }
425 3682 : plumed_assert( !getPntrToArgument(i)->hasDerivatives() );
426 : }
427 : }
428 : // The prefactor for average and sum is set here so the number of input scalars is guaranteed to be correct
429 794 : myfunc.setPrefactor( this, 1.0 );
430 2307 : return nelements;
431 : }
432 :
433 : template <class T>
434 11264 : void FunctionOfVector<T>::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
435 11264 : if( task_reducing_actions.size()==0 ) {
436 : return;
437 : }
438 1232 : if( !myfunc.allComponentsRequired( getArguments(), task_reducing_actions ) ) {
439 749 : task_reducing_actions.push_back(this);
440 : }
441 : }
442 :
443 : template <class T>
444 5541 : void FunctionOfVector<T>::runSingleTaskCalculation( const Value* arg, ActionWithValue* action, T& f ) {
445 : // This is used if we are doing sorting actions on a single vector
446 5541 : unsigned nv = arg->getNumberOfValues();
447 5541 : std::vector<double> args( nv );
448 8198467 : for(unsigned i=0; i<nv; ++i) {
449 8192926 : args[i] = arg->get(i);
450 : }
451 5541 : std::vector<double> vals( action->getNumberOfComponents() );
452 5541 : Matrix<double> derivatives( action->getNumberOfComponents(), nv );
453 5541 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(action);
454 5541 : plumed_assert( aa );
455 5541 : f.calc( aa, args, vals, derivatives );
456 11498 : for(unsigned i=0; i<vals.size(); ++i) {
457 5957 : action->copyOutput(i)->set( vals[i] );
458 : }
459 : // Return if we are not computing derivatives
460 5541 : if( action->doNotCalculateDerivatives() ) {
461 : return;
462 : }
463 : // Now set the derivatives
464 198059 : for(unsigned j=0; j<nv; ++j) {
465 388720 : for(unsigned i=0; i<vals.size(); ++i) {
466 195074 : action->copyOutput(i)->setDerivative( j, derivatives(i,j) );
467 : }
468 : }
469 : }
470 :
471 : template <class T>
472 185512 : void FunctionOfVector<T>::calculate() {
473 : // Everything is done elsewhere
474 185512 : if( actionInChain() ) {
475 : return;
476 : }
477 : // This is done if we are calculating a function of multiple cvs
478 80254 : if( !doAtEnd ) {
479 74713 : runAllTasks();
480 : }
481 : // This is used if we are doing sorting actions on a single vector
482 5541 : else if( !myfunc.doWithTasks() ) {
483 5541 : runSingleTaskCalculation( getPntrToArgument(0), this, myfunc );
484 : }
485 : }
486 :
487 : }
488 : }
489 : #endif
|
