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_FunctionOfMatrix_h
23 : #define __PLUMED_function_FunctionOfMatrix_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "FunctionSetup.h"
28 : #include "FunctionOfVector.h"
29 : #include "Custom.h"
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 : template <class CV, typename myPTM=defaultPTM>
35 : class FunctionOfMatrix : public ActionWithVector {
36 : public:
37 : using input_type = FunctionData<CV>;
38 : using mytype = FunctionOfMatrix<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 : /// Number of scalars that appear in the input
48 : std::size_t nscalars;
49 : /// Used to hold the list of tasks we are running
50 : std::vector<unsigned> active_tasks;
51 : /// Get the number of arguments the function uses
52 : unsigned getNumberOfFunctionArguments() const ;
53 : public:
54 : static void registerKeywords(Keywords&);
55 : explicit FunctionOfMatrix(const ActionOptions&);
56 : /// Get the label to write in the graph
57 : std::string writeInGraph() const override ;
58 : /// Get the number of derivatives for this action
59 : unsigned getNumberOfDerivatives() override ;
60 : /// Resize the matrices
61 : void prepare() override ;
62 : void calculate() override ;
63 : void getNumberOfTasks( unsigned& ntasks ) override ;
64 : std::vector<unsigned>& getListOfActiveTasks() override ;
65 : void getArgumentInputData( std::vector<double>& inputdata ) const override ;
66 : void getArgumentInputData( std::vector<float>& inputdata ) const override ;
67 : static void performTask( std::size_t task_index,
68 : const FunctionData<CV>& actiondata,
69 : ParallelActionsInput& input,
70 : ParallelActionsOutput& output );
71 : /// Add some forces
72 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
73 : /// Get the indices of the forces
74 : static int getNumberOfValuesPerTask( std::size_t task_index, const FunctionData<CV>& actiondata );
75 : static void getForceIndices( std::size_t task_index,
76 : std::size_t colno,
77 : std::size_t ntotal_force,
78 : const FunctionData<CV>& actiondata,
79 : const ParallelActionsInput& input,
80 : ForceIndexHolder force_indices );
81 : };
82 :
83 : template <class CV, typename myPTM>
84 913 : void FunctionOfMatrix<CV,myPTM>::registerKeywords(Keywords& keys ) {
85 913 : ActionWithVector::registerKeywords(keys);
86 913 : std::string name = keys.getDisplayName();
87 913 : std::size_t und=name.find("_MATRIX");
88 913 : keys.setDisplayName( name.substr(0,und) );
89 1826 : keys.addInputKeyword("compulsory","ARG","scalar/matrix","the labels of the scalar and matrices that on which the function is being calculated elementwise");
90 1826 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of the output matrix to compute");
91 913 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
92 913 : 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");
93 913 : CV::registerKeywords( keys );
94 913 : if( keys.getDisplayName()=="SUM" ) {
95 0 : keys.setValueDescription("scalar","the sum of all the elements in the input matrix");
96 913 : } else if( keys.getDisplayName()=="HIGHEST" ) {
97 0 : keys.setValueDescription("scalar","the largest element of the input matrix");
98 913 : } else if( keys.getDisplayName()=="LOWEST" ) {
99 0 : keys.setValueDescription("scalar","the smallest element in the input matrix");
100 1826 : } else if( keys.outputComponentExists(".#!value") ) {
101 3240 : keys.setValueDescription("matrix","the matrix obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input matrix");
102 : }
103 913 : PTM::registerKeywords( keys );
104 913 : }
105 :
106 : template <class CV, typename myPTM>
107 144864 : unsigned FunctionOfMatrix<CV,myPTM>::getNumberOfFunctionArguments() const {
108 144864 : unsigned nargs=getNumberOfArguments();
109 144864 : if( getNumberOfMasks()>0 ) {
110 946 : return nargs - getNumberOfMasks();
111 : }
112 : return nargs;
113 : }
114 :
115 : template <class CV, typename myPTM>
116 438 : FunctionOfMatrix<CV,myPTM>::FunctionOfMatrix(const ActionOptions&ao):
117 : Action(ao),
118 : ActionWithVector(ao),
119 438 : taskmanager(this),
120 438 : argstart(0),
121 438 : nscalars(0) {
122 : // Check if first argument is grid
123 438 : if( getPntrToArgument(0)->getRank()>0 && getPntrToArgument(0)->hasDerivatives() ) {
124 0 : argstart=1;
125 : }
126 438 : if( getNumberOfArguments()==argstart ) {
127 0 : error("no arguments specified");
128 : }
129 :
130 438 : if( getPntrToArgument(argstart)->getRank()!=2 ) {
131 0 : error("first argument to this action must be a matrix");
132 : }
133 :
134 : // Get the number of arguments
135 : unsigned nargs = getNumberOfArguments();
136 : int nmasks = getNumberOfMasks();
137 438 : if( nargs>=static_cast<unsigned>(nmasks) && nmasks>0 ) {
138 52 : nargs = nargs - nmasks;
139 : }
140 : // Get the shape of the output
141 440 : std::vector<std::size_t> shape( 2 );
142 438 : shape[0] = getPntrToArgument(argstart)->getShape()[0];
143 438 : shape[1] = getPntrToArgument(argstart)->getShape()[1];
144 958 : for(unsigned i=argstart+1; i<nargs; ++i) {
145 520 : if( getPntrToArgument(i)->getRank()==0 ) {
146 2 : nscalars++;
147 518 : } else if( getPntrToArgument(i)->getRank()==2 ) {
148 518 : if( getPntrToArgument(i)->getShape()[0]!=shape[0] || getPntrToArgument(i)->getShape()[1]!=shape[1] ) {
149 0 : error("mismatch between sizes of input arguments");
150 518 : } else if( nscalars>0 ) {
151 0 : error("scalars should be specified in argument list after all matrices");
152 : }
153 : } else {
154 0 : error("input arguments should be matrices or scalars");
155 : }
156 : }
157 438 : if( nmasks>0 ) {
158 52 : if( getPntrToArgument(getNumberOfArguments()-nmasks)->getShape()[0]!=shape[0] ||
159 52 : getPntrToArgument(getNumberOfArguments()-nmasks)->getShape()[1]!=shape[1] ) {
160 0 : error("input mask has wrong size");
161 : }
162 : }
163 :
164 : // Check if the output matrix is symmetric
165 438 : bool symmetric=true;
166 1448 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
167 1010 : if( getPntrToArgument(i)->getRank()==2 ) {
168 1008 : if( !getPntrToArgument(i)->isSymmetric() ) {
169 832 : symmetric=false;
170 : }
171 : }
172 : }
173 : // Setup the values
174 : // Get the function data from the parallel task manager, to avoid copies
175 438 : auto & myfunc = taskmanager.getActionInput();
176 438 : myfunc.argstart = argstart;
177 438 : myfunc.nscalars = nscalars;
178 438 : FunctionData<CV>::setup( myfunc.f, keywords.getOutputComponents(), shape, false, this );
179 : // Copy the fact that this is a symmetric matrix if the input matrices are all symmetric
180 1558 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
181 1122 : getPntrToComponent(i)->setSymmetric( symmetric );
182 : }
183 436 : taskmanager.setupParallelTaskManager( getNumberOfFunctionArguments() - argstart,
184 : nscalars );
185 444 : }
186 :
187 : template <class CV, typename myPTM>
188 0 : std::string FunctionOfMatrix<CV,myPTM>::writeInGraph() const {
189 0 : std::size_t und = getName().find_last_of("_");
190 0 : return getName().substr(0,und);
191 : }
192 :
193 : template <class CV, typename myPTM>
194 27635 : unsigned FunctionOfMatrix<CV,myPTM>::getNumberOfDerivatives() {
195 : unsigned nder=0;
196 136972 : for(unsigned i=argstart; i<getNumberOfFunctionArguments(); ++i) {
197 109337 : nder += getPntrToArgument(i)->getNumberOfStoredValues();
198 : }
199 27635 : return nder;
200 : }
201 :
202 : template <class CV, typename myPTM>
203 3730 : void FunctionOfMatrix<CV,myPTM>::prepare() {
204 3730 : std::vector<std::size_t> shape(getPntrToArgument(argstart)->getShape());
205 9809 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
206 6079 : Value* myval = getPntrToComponent(i);
207 6079 : if( myval->getRank()==2 && (myval->getShape()[0]!=shape[0] || myval->getShape()[1]!=shape[1]) ) {
208 18 : myval->setShape(shape);
209 : }
210 : }
211 3730 : ActionWithVector::prepare();
212 3730 : active_tasks.resize(0);
213 3730 : }
214 :
215 : template <class CV, typename myPTM>
216 4164 : void FunctionOfMatrix<CV,myPTM>::getNumberOfTasks( unsigned& ntasks ) {
217 4164 : ntasks=getPntrToComponent(0)->getNumberOfStoredValues();
218 4164 : }
219 :
220 : template <class CV, typename myPTM>
221 6199 : std::vector<unsigned>& FunctionOfMatrix<CV,myPTM>::getListOfActiveTasks() {
222 6199 : if( active_tasks.size()>0 ) {
223 2476 : return active_tasks;
224 : }
225 :
226 : Value* myarg = NULL;
227 3723 : if( getNumberOfMasks()>0 ) {
228 343 : myarg = getPntrToArgument(getNumberOfArguments()-getNumberOfMasks());
229 : } else {
230 3380 : myarg = getPntrToArgument(argstart);
231 : }
232 : unsigned atsize = 0;
233 3723 : unsigned nrows = myarg->getShape()[0];
234 1163596 : for(unsigned i=0; i<nrows; ++i) {
235 2319746 : atsize += myarg->getRowLength(i);
236 : }
237 3723 : active_tasks.resize( atsize );
238 :
239 1163596 : for(unsigned i=0, base=0,k=0; i<nrows; ++i) {
240 1159873 : unsigned ncols = myarg->getRowLength(i);
241 34452216 : for(unsigned j=0; j<ncols; ++j) {
242 33292343 : active_tasks[k] = base+j;
243 33292343 : ++k;
244 : }
245 1159873 : base += myarg->getNumberOfColumns();
246 : }
247 3723 : if( getNumberOfMasks()>0 && doNotCalculateDerivatives() ) {
248 : return active_tasks;
249 : }
250 : // All the matrices input have to have the same sparsity pattern.
251 : // I can do everything I want to do with this limitation. If
252 : // anyone wants to make this smarter in the future they can
253 : // This check was guarded by a misspelled "#ifndef DNDEBUG", so it has always been active.
254 : // It is kept active also in optimized builds: it costs ~1% of a CONTACT_MATRIX->CUSTOM pipeline
255 : // and it catches inputs that would otherwise give wrong results silently.
256 10148 : for(unsigned k=argstart; k<getNumberOfArguments(); ++k) {
257 6612 : if( getPntrToArgument(k)->getRank()!=2 ) {
258 592 : continue ;
259 : }
260 6020 : if( getNumberOfMasks()>0 && getPntrToArgument(k)->isConstant() ) {
261 10 : continue ;
262 : }
263 2931255 : for(unsigned i=0; i<nrows; ++i) {
264 2925245 : unsigned ncols = myarg->getRowLength(i);
265 2925245 : if( getNumberOfMasks()>0 && ncols==0 ) {
266 522 : continue;
267 : }
268 2924723 : plumed_massert( ncols==getPntrToArgument(k)->getRowLength(i), "failing in " + getLabel() );
269 65037366 : for(unsigned j=0; j<ncols; ++j) {
270 62112643 : plumed_assert( myarg->getRowIndex(i,j)==getPntrToArgument(k)->getRowIndex(i,j) );
271 : }
272 : }
273 : }
274 : return active_tasks;
275 : }
276 :
277 : template <class CV, typename myPTM>
278 3728 : void FunctionOfMatrix<CV,myPTM>::getArgumentInputData( std::vector<double>& inputdata ) const {
279 : int nmasks = getNumberOfMasks();
280 3728 : unsigned nargs = getNumberOfFunctionArguments();
281 :
282 3728 : const Value* myval = getConstPntrToComponent(0);
283 3728 : std::size_t ntasks = myval->getNumberOfStoredValues();
284 3728 : std::size_t ndata = static_cast<std::size_t>(nargs-argstart)*ntasks;
285 3728 : if( inputdata.size()!=ndata ) {
286 578 : inputdata.resize( ndata );
287 : }
288 : const Value* maskarg=NULL;
289 3728 : if( getNumberOfMasks()>0 ) {
290 343 : maskarg = getPntrToArgument(getNumberOfArguments()-getNumberOfMasks());
291 : }
292 :
293 10464 : for(unsigned j=argstart; j<nargs; ++j) {
294 : const Value* jarg = getPntrToArgument(j);
295 6736 : if( jarg->getRank()==0 ) {
296 592 : double val = jarg->get();
297 2403 : for(unsigned i=0; i<myval->getShape()[0]; ++i) {
298 1811 : unsigned colbase=i*myval->getNumberOfColumns();
299 9355 : for(unsigned k=0; k<myval->getRowLength(i); ++k) {
300 5733 : inputdata[(nargs-argstart)*(colbase+k) + j-argstart] = val;
301 : }
302 : }
303 6144 : } else if( nmasks>0 && jarg!=maskarg ) {
304 56565 : for(unsigned i=0; i<myval->getShape()[0]; ++i) {
305 56091 : unsigned jcolbase = i*jarg->getShape()[1];
306 56091 : unsigned vcolbase = i*myval->getNumberOfColumns();
307 16206256 : for(unsigned k=0; k<myval->getRowLength(i); ++k) {
308 8047037 : inputdata[(nargs-argstart)*(vcolbase+k) + j-argstart] = jarg->get(jcolbase+myval->getRowIndex(i,k),true);
309 : }
310 : }
311 : } else {
312 2933790 : for(unsigned i=0; i<jarg->getShape()[0]; ++i) {
313 2928120 : unsigned colbase=i*jarg->getNumberOfColumns();
314 128882378 : for(unsigned k=0; k<jarg->getRowLength(i); ++k) {
315 61513069 : inputdata[(nargs-argstart)*(colbase+k) + j-argstart] = jarg->get(colbase+k,false);
316 : }
317 : }
318 : }
319 : }
320 3728 : }
321 :
322 : template <class CV, typename myPTM>
323 0 : void FunctionOfMatrix<CV,myPTM>::getArgumentInputData( std::vector<float>& inputdata ) const {
324 : int nmasks = getNumberOfMasks();
325 0 : unsigned nargs = getNumberOfFunctionArguments();
326 :
327 0 : const Value* myval = getConstPntrToComponent(0);
328 0 : std::size_t ntasks = myval->getNumberOfStoredValues();
329 0 : std::size_t ndata = static_cast<std::size_t>(nargs-argstart)*ntasks;
330 0 : if( inputdata.size()!=ndata ) {
331 0 : inputdata.resize( ndata );
332 : }
333 : const Value* maskarg=NULL;
334 0 : if( getNumberOfMasks()>0 ) {
335 0 : maskarg = getPntrToArgument(getNumberOfArguments()-getNumberOfMasks());
336 : }
337 :
338 0 : for(unsigned j=argstart; j<nargs; ++j) {
339 : const Value* jarg = getPntrToArgument(j);
340 0 : if( jarg->getRank()==0 ) {
341 0 : double val = jarg->get();
342 0 : for(unsigned i=0; i<myval->getShape()[0]; ++i) {
343 0 : unsigned colbase=i*myval->getNumberOfColumns();
344 0 : for(unsigned k=0; k<myval->getRowLength(i); ++k) {
345 0 : inputdata[(nargs-argstart)*(colbase+k) + j-argstart] = val;
346 : }
347 : }
348 0 : } else if( nmasks>0 && jarg!=maskarg ) {
349 0 : for(unsigned i=0; i<myval->getShape()[0]; ++i) {
350 0 : unsigned jcolbase = i*jarg->getShape()[1];
351 0 : unsigned vcolbase = i*myval->getNumberOfColumns();
352 0 : for(unsigned k=0; k<myval->getRowLength(i); ++k) {
353 0 : inputdata[(nargs-argstart)*(vcolbase+k) + j-argstart] = jarg->get(jcolbase+myval->getRowIndex(i,k),true);
354 : }
355 : }
356 : } else {
357 0 : for(unsigned i=0; i<jarg->getShape()[0]; ++i) {
358 0 : unsigned colbase=i*jarg->getNumberOfColumns();
359 0 : for(unsigned k=0; k<jarg->getRowLength(i); ++k) {
360 0 : inputdata[(nargs-argstart)*(colbase+k) + j-argstart] = jarg->get(colbase+k,false);
361 : }
362 : }
363 : }
364 : }
365 0 : }
366 :
367 : template <class CV, typename myPTM>
368 3728 : void FunctionOfMatrix<CV,myPTM>::calculate() {
369 : Value* myarg = NULL;
370 3728 : if( getNumberOfMasks()>0 ) {
371 343 : myarg = getPntrToArgument(getNumberOfArguments()-getNumberOfMasks());
372 : } else {
373 3385 : myarg = getPntrToArgument(argstart);
374 : }
375 : // Copy bookeeping arrays from input matrices to output matrices
376 9805 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
377 6077 : getPntrToComponent(i)->copyBookeepingArrayFromArgument( myarg );
378 : }
379 3728 : taskmanager.setupParallelTaskManager( getNumberOfFunctionArguments()-argstart, nscalars );
380 3728 : taskmanager.runAllTasks();
381 3728 : }
382 :
383 : template <class CV, typename myPTM>
384 32516764 : void FunctionOfMatrix<CV,myPTM>::performTask( std::size_t task_index,
385 : const FunctionData<CV>& actiondata,
386 : ParallelActionsInput& input,
387 : ParallelActionsOutput& output ) {
388 32516764 : auto funcout = FunctionOutput::create( input.ncomponents,
389 : output.values.data(),
390 32516764 : input.nderivatives_per_scalar,
391 : output.derivatives.data() );
392 35283300 : CV::calc( actiondata.f,
393 32516764 : input.noderiv,
394 32516764 : View<const double>( input.inputdata + task_index*input.nderivatives_per_scalar,
395 : input.nderivatives_per_scalar ),
396 : funcout );
397 32516764 : }
398 :
399 : template <class CV, typename myPTM>
400 2471 : void FunctionOfMatrix<CV,myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
401 2471 : taskmanager.applyForces( outforces );
402 2471 : }
403 :
404 : template <class CV, typename myPTM>
405 3934037 : int FunctionOfMatrix<CV,myPTM>::getNumberOfValuesPerTask( std::size_t task_index, const FunctionData<CV>& actiondata ) {
406 3934037 : return 1;
407 : }
408 :
409 : template <class CV, typename myPTM>
410 3837986 : void FunctionOfMatrix<CV,myPTM>::getForceIndices( std::size_t task_index,
411 : std::size_t colno,
412 : std::size_t ntotal_force,
413 : const FunctionData<CV>& actiondata,
414 : const ParallelActionsInput& input,
415 : ForceIndexHolder force_indices ) {
416 : // The force indices are found in the same way as in FunctionOfVector so we reuse that function here
417 3837986 : FunctionOfVector<CV>::getForceIndices( task_index, colno, ntotal_force, actiondata, input, force_indices );
418 3837986 : }
419 :
420 :
421 :
422 : /// Explicit specialization defined in the corresponding .cpp file; declared
423 : /// here so that every translation unit instantiating the template sees it.
424 : template <> std::string FunctionOfMatrix<Custom>::writeInGraph() const;
425 : }
426 : }
427 : #endif
|
