Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2017 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_gridtools_KDE_h
23 : #define __PLUMED_gridtools_KDE_h
24 :
25 : #include "ActionWithGrid.h"
26 : #include "SumOfKernels.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 : #include "core/PbcAction.h"
30 : #include "tools/HistogramBead.h"
31 : #include "core/ParallelTaskManager.h"
32 : #include "tools/SwitchingFunction.h"
33 : #include "tools/Matrix.h"
34 :
35 : namespace PLMD {
36 : namespace gridtools {
37 :
38 : template <class K, class P, class G>
39 : class KDEHelper {
40 : public:
41 : G g;
42 : bool fixed_width;
43 : std::size_t maxkernels;
44 : SumOfKernels<K, P> kernelsum;
45 : std::vector<unsigned> nneigh;
46 : std::vector<std::size_t> nkernels_per_point;
47 : std::vector<std::size_t> kernels_for_gridpoint;
48 : static void registerKeywords( Keywords& keys );
49 : static void read( KDEHelper<K,P,G>& func,
50 : ActionWithArguments* action,
51 : const std::vector<Value*>& args,
52 : GridCoordinatesObject& gridobject,
53 : std::vector<std::size_t>& shape,
54 : function::FunctionOptions& options );
55 : static void readKernelParameters( std::string& value, ActionWithArguments* action, const std::string& outlab, bool rerequestargs );
56 : static void addArgument( const std::string& value, ActionWithArguments* action );
57 : static void setupGridBounds( KDEHelper<K,P,G>& func, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval );
58 : static void transferParamsToKernel( const std::vector<double>& argval, KDEHelper<K,P,G>& func, GridCoordinatesObject& gridobject, bool updateNeighborsOnEachKernel, std::size_t nkernels, unsigned kval, K& kp );
59 : static void transferKernels( KDEHelper<K,P,G>& func, const std::vector<Value*>& args, GridCoordinatesObject& gridobject );
60 : };
61 :
62 : template <class K, class P, class G>
63 323 : void KDEHelper<K,P,G>::registerKeywords( Keywords& keys ) {
64 282 : SumOfKernels<K,P>::registerKeywords( keys );
65 323 : G::registerKeywords( keys );
66 323 : }
67 :
68 : template <class K, class P, class G>
69 89 : void KDEHelper<K,P,G>::read( KDEHelper<K,P,G>& func,
70 : ActionWithArguments* action,
71 : const std::vector<Value*>& args,
72 : GridCoordinatesObject& gridobject,
73 : std::vector<std::size_t>& shape,
74 : function::FunctionOptions& options ) {
75 89 : func.nneigh.resize( args.size() );
76 66 : SumOfKernels<K,P>::read( func.kernelsum, action, args, options );
77 89 : G::readBandwidthAndHeight( func.kernelsum.params, action );
78 339 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
79 250 : if( (action->getPntrToArgument(0))->getNumberOfValues()==0 && (action->getPntrToArgument(i))->isConstant() ) {
80 1 : continue;
81 : }
82 :
83 249 : if( (action->getPntrToArgument(0))->getNumberOfValues()!=1 || (action->getPntrToArgument(0))->getNumberOfValues()!=1 ) {
84 173 : if( (action->getPntrToArgument(0))->getRank()!=(action->getPntrToArgument(i))->getRank() ) {
85 0 : action->error("mismatch between ranks of input actions");
86 : }
87 373 : for(unsigned j=0; j<(action->getPntrToArgument(0))->getRank(); ++j) {
88 200 : if( (action->getPntrToArgument(0))->getShape()[j]!=(action->getPntrToArgument(i))->getShape()[j] ) {
89 0 : action->error("mismatch between shapes of input actions");
90 : }
91 : }
92 : }
93 : }
94 89 : G::readGridParameters( func.g, action, gridobject, shape );
95 89 : }
96 :
97 : template <class K, class P, class G>
98 86 : void KDEHelper<K,P,G>::setupGridBounds( KDEHelper<K,P,G>& func, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval ) {
99 : // Setup the grid boundaries on first step
100 86 : G::setupGridBounds( func.g, box, gridobject, args, myval );
101 : // Check if the bandwidth changes during the simulation
102 86 : func.fixed_width = false;
103 86 : if( K::bandwidthIsConstant( gridobject.getDimension(), args ) && K::bandwidthsAllSame( gridobject.getDimension(), args ) ) {
104 0 : K myk;
105 86 : std::vector<double> myargs( args.size() );
106 409 : for(unsigned j=0; j<args.size(); ++j) {
107 323 : myargs[j] = args[j]->get(0);
108 : }
109 86 : K::setKernelAndCheckHeight( myk, gridobject.getDimension(), myargs );
110 86 : G::getDiscreteSupport( func.g, func.kernelsum.params, myk, func.nneigh, gridobject );
111 86 : func.fixed_width = true;
112 0 : }
113 172 : if( gridobject.getGridType()=="fibonacci" ) {
114 : return;
115 : }
116 : // Set the periodicity of the parameters
117 181 : for(unsigned i=0; i<gridobject.getDimension(); ++i) {
118 192 : P::setArgumentDomain( i, func.kernelsum.params, gridobject.getGridSpacing()[i], gridobject.isPeriodic(i), gridobject.getMin()[i], gridobject.getMax()[i] );
119 : }
120 : }
121 :
122 : template <class K, class P, class G>
123 149 : void KDEHelper<K,P,G>::readKernelParameters( std::string& value, ActionWithArguments* action, const std::string& outlab, bool rerequestargs ) {
124 : std::size_t dot = value.find_first_of('.');
125 149 : ActionWithValue* av = action->plumed.getActionSet().selectWithLabel<ActionWithValue*>( value.substr(0,dot) );
126 149 : if( !av ) {
127 108 : std::string matstr, vals = "VALUES=" + value;
128 108 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
129 : std::string nr, nc;
130 8 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
131 8 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
132 16 : matstr = " NROWS=" + nr + " NCOLS=" + nc;
133 : }
134 788438 : for(unsigned i=1; i<(action->getPntrToArgument(0))->getNumberOfValues(); ++i) {
135 1576660 : vals += "," + value;
136 : }
137 108 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + ": CONSTANT " + vals + matstr ), false );
138 216 : value = action->getLabel() + outlab;
139 : } else {
140 : Value* myval;
141 41 : if( dot!=std::string::npos ) {
142 1 : myval = av->copyOutput( value );
143 : } else {
144 40 : if( av->getNumberOfComponents()>1 ) {
145 0 : action->error("problem reading argument " + value );
146 : }
147 40 : myval = av->copyOutput(0);
148 : }
149 41 : if( myval->getRank()==0 ) {
150 : std::string nvals;
151 0 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
152 : std::string nr, nc;
153 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
154 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
155 0 : nvals = nr + "," + nc;
156 : } else {
157 0 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), nvals );
158 : }
159 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + "_ones: ONES SIZE=" + nvals ), false );
160 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + ": CUSTOM ARG=" + action->getLabel() + outlab + "_ones," + value ), false );
161 0 : value = action->getLabel() + outlab;
162 : }
163 : }
164 149 : if( !rerequestargs ) {
165 : return;
166 : }
167 125 : KDEHelper<K,P,G>::addArgument( value, action );
168 : }
169 :
170 : template <class K, class P, class G>
171 197 : void KDEHelper<K,P,G>::addArgument( const std::string& value, ActionWithArguments* action ) {
172 197 : std::size_t dot = value.find_first_of(".");
173 197 : ActionWithValue* av = action->plumed.getActionSet().selectWithLabel<ActionWithValue*>( value.substr(0,dot) );
174 197 : plumed_assert( av );
175 197 : std::vector<Value*> args( action->getArguments() );
176 197 : args.push_back( av->copyOutput(0) );
177 197 : action->requestArguments( args );
178 197 : }
179 :
180 : template <class K, class P, class G>
181 624237 : void KDEHelper<K,P,G>::transferParamsToKernel( const std::vector<double>& argval, KDEHelper<K,P,G>& func, GridCoordinatesObject& gridobject, bool updateNeighborsOnEachKernel, std::size_t nkernels, unsigned kval, K& kp ) {
182 : // This sets the kernel parameters for the Kth kernel and checks that we want
183 : // to consider it
184 624237 : if( !K::setKernelAndCheckHeight( kp, gridobject.getDimension(), argval ) ) {
185 194981 : return;
186 : }
187 : // If the widths of each kernel are not all the same then get the discrete support
188 429256 : if( updateNeighborsOnEachKernel ) {
189 0 : G::getDiscreteSupport( func.g, func.kernelsum.params, kp, func.nneigh, gridobject );
190 : }
191 :
192 : // Now get the grid points for this particular kernel
193 : unsigned num_neigh;
194 : std::vector<unsigned> neighbors;
195 429256 : G::getNeighbors( func.kernelsum.params, kp, gridobject, func.nneigh, num_neigh, neighbors );
196 :
197 : // And transfer the neighbor information to the holders
198 60663074 : for(unsigned j=0; j<num_neigh; ++j) {
199 60233818 : func.kernels_for_gridpoint[ neighbors[j]*nkernels + func.nkernels_per_point[neighbors[j]] ] = kval;
200 60233818 : func.nkernels_per_point[ neighbors[j] ]++;
201 : }
202 : }
203 :
204 : template <class K, class P, class G>
205 2344 : void KDEHelper<K,P,G>::transferKernels( KDEHelper<K,P,G>& func, const std::vector<Value*>& args, GridCoordinatesObject& gridobject ) {
206 : // Resize the kernel sum if we need to
207 : // Number of kernels is determined based on sparsity pattern of matrix input as matrix of heights
208 2344 : std::size_t nkernels = args[args.size()-1]->getNumberOfStoredValues();
209 2344 : if( func.kernelsum.kernelParams.size()!=nkernels ) {
210 90 : func.kernelsum.kernelParams.resize( nkernels );
211 : }
212 : // And resize the grid counters if we need to
213 2344 : std::size_t ngp = gridobject.getNumberOfPoints();
214 2344 : if( func.nkernels_per_point.size()!=ngp ) {
215 86 : func.nkernels_per_point.resize( ngp );
216 86 : func.kernels_for_gridpoint.resize( ngp*nkernels );
217 : }
218 : std::fill( func.nkernels_per_point.begin(), func.nkernels_per_point.end(), 0 );
219 :
220 2344 : bool updateNeighborsOnEachKernel = !func.fixed_width;
221 2344 : if( !func.fixed_width && K::bandwidthsAllSame( gridobject.getDimension(), args ) ) {
222 0 : G::getDiscreteSupport( func.g, func.kernelsum.params, func.kernelsum.kernelParams[0], func.nneigh, gridobject );
223 : updateNeighborsOnEachKernel = false;
224 : }
225 :
226 2344 : std::vector<double> argval( args.size() );
227 2344 : if( args[args.size()-1]->getRank()==2 ) {
228 33 : const unsigned nc = args[args.size()-1]->getShape()[1];
229 33 : const unsigned nrows = args[args.size()-1]->getShape()[0];
230 33 : const unsigned ncs = args[args.size()-1]->getNumberOfColumns();
231 5577 : for(unsigned i=0; i<nrows; ++i) {
232 5544 : unsigned ncols = args[args.size()-1]->getRowLength(i);
233 28486 : for(unsigned k=0; k<args.size(); ++k) {
234 40240 : plumed_massert( args[k]->isConstant() || ncols==args[k]->getRowLength(i), "all input matrices must have same sparsity pattern" );
235 : }
236 466519 : for(unsigned j=0; j<ncols; ++j) {
237 460975 : unsigned jind = args[args.size()-1]->getRowIndex( i, j );
238 2233088 : for(unsigned k=0; k<args.size(); ++k) {
239 1772113 : if( jind==args[k]->getRowIndex( i, j ) ) {
240 1331589 : argval[k] = args[k]->get( i*ncs+ j, false );
241 : } else {
242 440524 : argval[k] = args[k]->get( i*nc + jind );
243 : }
244 : }
245 460975 : KDEHelper<K,P,G>::transferParamsToKernel( argval, func, gridobject, updateNeighborsOnEachKernel, nkernels, i*ncs+j, func.kernelsum.kernelParams[i*ncs+j] );
246 : }
247 : }
248 : } else {
249 165573 : for(unsigned i=0; i<nkernels; ++i) {
250 : // Transfer the kernel parameters to local vector of doubles
251 871820 : for(unsigned j=0; j<args.size(); ++j) {
252 708558 : argval[j] = args[j]->get(i,false);
253 : }
254 163262 : KDEHelper<K,P,G>::transferParamsToKernel( argval, func, gridobject, updateNeighborsOnEachKernel, nkernels, i, func.kernelsum.kernelParams[i] );
255 : }
256 : }
257 : // Get the maximum number of kernels for any given grid point (used for resizing derivatives)
258 2344 : func.maxkernels = 0;
259 577517 : for(unsigned i=0; i<ngp; ++i) {
260 575173 : if( func.nkernels_per_point[i]>func.maxkernels ) {
261 3398 : func.maxkernels = func.nkernels_per_point[i];
262 : }
263 : }
264 2344 : }
265 :
266 : template <class K, class P, class G>
267 : class KDE : public ActionWithGrid {
268 : public:
269 : using input_type = KDEHelper<K, P, G>;
270 : using PTM = ParallelTaskManager<KDE<K,P,G>>;
271 : private:
272 : bool firststep;
273 : /// The parallel task manager
274 : PTM taskmanager;
275 : GridCoordinatesObject gridobject;
276 : public:
277 : static void registerKeywords( Keywords& keys );
278 : explicit KDE(const ActionOptions&ao);
279 : std::vector<std::string> getGridCoordinateNames() const override ;
280 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
281 : unsigned getNumberOfDerivatives() override;
282 : int checkTaskIsActive( const unsigned& itask ) const override ;
283 : void prepare() override ;
284 : void calculate() override ;
285 : void getInputData( std::vector<double>& inputdata ) const override ;
286 : static void performTask( std::size_t task_index,
287 : const KDEHelper<K, P, G>& actiondata,
288 : ParallelActionsInput& input,
289 : ParallelActionsOutput& output );
290 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
291 : static int getNumberOfValuesPerTask( std::size_t task_index,
292 : const KDEHelper<K, P, G>& actiondata );
293 : static void getForceIndices( std::size_t task_index,
294 : std::size_t colno,
295 : std::size_t ntotal_force,
296 : const KDEHelper<K, P, G>& actiondata,
297 : const ParallelActionsInput& input,
298 : ForceIndexHolder force_indices );
299 : };
300 :
301 : template <class K, class P, class G>
302 323 : void KDE<K,P,G>::registerKeywords( Keywords& keys ) {
303 323 : ActionWithGrid::registerKeywords( keys );
304 646 : keys.addInputKeyword("compulsory","ARG","scalar/vector/matrix","the label for the value that should be used to construct the histogram");
305 323 : KDEHelper<K,P,G>::registerKeywords( keys );
306 : // Keywords for spherical KDE
307 323 : keys.add("hidden","MASKED_INPUT_ALLOWED","turns on that you are allowed to use masked inputs ");
308 646 : keys.setValueDescription("grid","a function on a grid that was obtained by doing a Kernel Density Estimation using the input arguments");
309 646 : if( keys.getDisplayName()!="SPHERICAL_KDE" ) {
310 622 : keys.setDisplayName("KDE");
311 : }
312 323 : PTM::registerKeywords( keys );
313 323 : }
314 :
315 : template <class K, class P, class G>
316 89 : KDE<K,P,G>::KDE(const ActionOptions&ao):
317 : Action(ao),
318 : ActionWithGrid(ao),
319 89 : firststep(true),
320 89 : taskmanager(this) {
321 :
322 89 : std::vector<std::size_t> shape( getNumberOfArguments() );
323 89 : unsigned numberOfKernels=getPntrToArgument(0)->getNumberOfValues();
324 142 : for(unsigned i=1; i<shape.size(); ++i) {
325 53 : if( numberOfKernels!=getPntrToArgument(i)->getNumberOfValues() ) {
326 0 : error("mismatch between numbers of values in input arguments");
327 : }
328 : }
329 :
330 : function::FunctionOptions foptions;
331 89 : KDEHelper<K,P,G>::read( taskmanager.getActionInput(), this, getArguments(), gridobject, shape, foptions );
332 89 : addValueWithDerivatives( shape );
333 89 : setNotPeriodic();
334 89 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
335 89 : }
336 :
337 : template <class K, class P, class G>
338 7250 : unsigned KDE<K,P,G>::getNumberOfDerivatives() {
339 7250 : return gridobject.getDimension();
340 : }
341 :
342 : template <class K, class P, class G>
343 141 : std::vector<std::string> KDE<K,P,G>::getGridCoordinateNames() const {
344 141 : std::vector<std::string> names( gridobject.getDimension() );
345 354 : for(unsigned i=0; i<names.size(); ++i) {
346 : names[i] = getPntrToArgument(i)->getName();
347 : }
348 141 : return names;
349 0 : }
350 :
351 : template <class K, class P, class G>
352 279165 : const GridCoordinatesObject& KDE<K,P,G>::getGridCoordinatesObject() const {
353 279165 : return gridobject;
354 : }
355 :
356 : template <class K, class P, class G>
357 575173 : int KDE<K,P,G>::checkTaskIsActive( const unsigned& itask ) const {
358 575173 : if( taskmanager.getActionInput().nkernels_per_point[itask]>0 ) {
359 177913 : return 1;
360 : }
361 : return -1;
362 : }
363 :
364 : template <class K, class P, class G>
365 2344 : void KDE<K,P,G>::prepare() {
366 2344 : ActionWithVector::prepare();
367 2344 : std::size_t nkernels = getPntrToArgument(0)->getNumberOfValues();
368 7297 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
369 : Value* myarg = getPntrToArgument(i);
370 4953 : if( myarg->getNumberOfValues()!=nkernels ) {
371 1 : if( myarg->isConstant() && myarg->getNumberOfValues()==1 ) {
372 1 : myarg->reshapeConstantValue( getPntrToArgument(0)->getShape() );
373 : } else {
374 0 : plumed_merror("found mismatched numbers of arguments in input");
375 : }
376 : }
377 : }
378 2344 : }
379 :
380 : template <class K, class P, class G>
381 2344 : void KDE<K,P,G>::calculate() {
382 2344 : if( firststep ) {
383 86 : PbcAction* bv = plumed.getActionSet().template selectWithLabel<PbcAction*>("Box");
384 86 : KDEHelper<K,P,G>::setupGridBounds( taskmanager.getActionInput(), bv->getPbc().getBox(), gridobject, getArguments(), getPntrToComponent(0) );
385 86 : firststep=false;
386 : }
387 2344 : KDEHelper<K,P,G>::transferKernels( taskmanager.getActionInput(), getArguments(), gridobject );
388 2344 : taskmanager.setupParallelTaskManager( getNumberOfArguments()*taskmanager.getActionInput().maxkernels, getNumberOfForceDerivatives() );
389 2344 : taskmanager.runAllTasks();
390 2344 : }
391 :
392 : template <class K, class P, class G>
393 2344 : void KDE<K,P,G>::getInputData( std::vector<double>& inputdata ) const {
394 2344 : std::size_t ndim = gridobject.getDimension();
395 2344 : std::size_t nstored = getConstPntrToComponent(0)->getNumberOfStoredValues();
396 2344 : std::vector<double> pos( ndim );
397 2344 : if( inputdata.size()!=nstored*ndim ) {
398 86 : inputdata.resize( ndim*nstored );
399 : }
400 :
401 577517 : for(unsigned i=0; i<nstored; ++i) {
402 575173 : gridobject.getGridPointCoordinates( i, pos );
403 1918145 : for(unsigned j=0; j<ndim; ++j) {
404 1342972 : inputdata[ i*ndim + j ] = pos[j];
405 : }
406 : }
407 2344 : }
408 :
409 : template <class K, class P, class G>
410 183915 : void KDE<K,P,G>::performTask( std::size_t task_index,
411 : const KDEHelper<K, P, G>& actiondata,
412 : ParallelActionsInput& input,
413 : ParallelActionsOutput& output ) {
414 : std::size_t ndim = actiondata.nneigh.size();
415 367830 : SumOfKernels<K,P>::calc( View<const std::size_t>( actiondata.kernels_for_gridpoint.data() + task_index*input.argstarts[1], actiondata.nkernels_per_point[task_index] ),
416 183915 : actiondata.kernelsum,
417 183915 : View<const double>(input.inputdata + task_index*actiondata.nneigh.size(),ndim),
418 : View<double>(output.values.data(), 1),
419 : View<double>(output.values.data()+1, ndim),
420 183915 : View<double>(output.derivatives.data(),actiondata.maxkernels*input.nargs) );
421 :
422 183915 : }
423 :
424 : template <class K, class P, class G>
425 612 : void KDE<K,P,G>::applyNonZeroRankForces( std::vector<double>& outforces ) {
426 612 : taskmanager.applyForces( outforces );
427 612 : }
428 :
429 : template <class K, class P, class G>
430 6002 : int KDE<K,P,G>::getNumberOfValuesPerTask( std::size_t task_index,
431 : const KDEHelper<K, P, G>& actiondata ) {
432 6002 : return 1;
433 : }
434 :
435 : template <class K, class P, class G>
436 6002 : void KDE<K,P,G>::getForceIndices( std::size_t task_index,
437 : std::size_t colno,
438 : std::size_t ntotal_force,
439 : const KDEHelper<K, P, G>& actiondata,
440 : const ParallelActionsInput& input,
441 : ForceIndexHolder force_indices ) {
442 6002 : force_indices.threadsafe_derivatives_end[0] = 0;
443 6002 : std::size_t nparams = K::getNumberOfParameters( actiondata.kernelsum.kernelParams[0] );
444 6002 : View<const std::size_t> kernellist( actiondata.kernels_for_gridpoint.data() + task_index*input.argstarts[1], actiondata.nkernels_per_point[task_index] );
445 4755078 : for(unsigned i=0; i<kernellist.size(); ++i) {
446 27660156 : for(unsigned j=0; j<nparams; ++j) {
447 22911080 : force_indices.indices[0][i*nparams+j] = input.argstarts[j] + kernellist[i];
448 : }
449 : }
450 6002 : force_indices.tot_indices[0] = kernellist.size()*nparams;
451 6002 : }
452 :
453 : }
454 : }
455 : #endif
|
