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 13 : class KDEHelper {
40 : public:
41 : G g;
42 : bool fixed_width;
43 : std::size_t maxkernels;
44 : std::size_t ndim{0};
45 : SumOfKernels<K, P> kernelsum;
46 : std::vector<unsigned> nneigh;
47 : std::vector<std::size_t> nkernels_per_point_v;
48 : std::size_t *nkernels_per_point{nullptr};
49 : std::vector<std::size_t> kernels_for_gridpoint_v;
50 : std::size_t *kernels_for_gridpoint{nullptr};
51 : static void registerKeywords( Keywords& keys );
52 : static void read( KDEHelper<K,P,G>& func,
53 : ActionWithArguments* action,
54 : const std::vector<Value*>& args,
55 : GridCoordinatesObject& gridobject,
56 : std::vector<std::size_t>& shape,
57 : function::FunctionOptions& options );
58 : static void readKernelParameters( std::string& value, ActionWithArguments* action, const std::string& outlab, bool rerequestargs );
59 : static void addArgument( const std::string& value, ActionWithArguments* action );
60 : static void setupGridBounds( KDEHelper<K,P,G>& func, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval );
61 : static void transferParamsToKernel( const std::vector<double>& argval, KDEHelper<K,P,G>& func, GridCoordinatesObject& gridobject, bool updateNeighborsOnEachKernel, std::size_t nkernels, unsigned kval, View<double>& params );
62 : static void transferKernels( KDEHelper<K,P,G>& func, const std::vector<Value*>& args, GridCoordinatesObject& gridobject );
63 : void update() {
64 2599 : ndim = nneigh.size();
65 2599 : nkernels_per_point = nkernels_per_point_v.data();
66 2599 : kernels_for_gridpoint = kernels_for_gridpoint_v.data();
67 : kernelsum.update();
68 : }
69 : void toACCDevice() const {
70 : #pragma acc enter data copyin(this[0:1],fixed_width,maxkernels,ndim, \
71 : nkernels_per_point[0:nkernels_per_point_v.size()], \
72 : kernels_for_gridpoint[0:kernels_for_gridpoint_v.size()])
73 : kernelsum.toACCDevice();
74 : }
75 : void removeFromACCDevice() const {
76 : kernelsum.removeFromACCDevice();
77 : #pragma acc exit data delete(kernels_for_gridpoint[0:kernels_for_gridpoint_v.size()], \
78 : nkernels_per_point[0:nkernels_per_point_v.size()], \
79 : ndim,maxkernels,fixed_width,this[0:1])
80 : }
81 : };
82 :
83 : template <class K, class P, class G>
84 443 : void KDEHelper<K,P,G>::registerKeywords( Keywords& keys ) {
85 362 : SumOfKernels<K,P>::registerKeywords( keys );
86 443 : G::registerKeywords( keys );
87 443 : }
88 :
89 : template <class K, class P, class G>
90 112 : void KDEHelper<K,P,G>::read( KDEHelper<K,P,G>& func,
91 : ActionWithArguments* action,
92 : const std::vector<Value*>& args,
93 : GridCoordinatesObject& gridobject,
94 : std::vector<std::size_t>& shape,
95 : function::FunctionOptions& options ) {
96 112 : func.nneigh.resize( args.size() );
97 84 : SumOfKernels<K,P>::read( func.kernelsum, action, args, options );
98 112 : G::readBandwidthAndHeight( func.kernelsum.params, action );
99 413 : for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
100 302 : if( (action->getPntrToArgument(0))->getNumberOfValues()==0 && (action->getPntrToArgument(i))->isConstant() ) {
101 1 : continue;
102 : }
103 :
104 301 : if( (action->getPntrToArgument(0))->getNumberOfValues()!=1 || (action->getPntrToArgument(0))->getNumberOfValues()!=1 ) {
105 201 : if( (action->getPntrToArgument(0))->getRank()!=(action->getPntrToArgument(i))->getRank() ) {
106 0 : action->error("mismatch between ranks of input actions");
107 : }
108 443 : for(unsigned j=0; j<(action->getPntrToArgument(0))->getRank(); ++j) {
109 242 : if( (action->getPntrToArgument(0))->getShape()[j]!=(action->getPntrToArgument(i))->getShape()[j] ) {
110 0 : action->error("mismatch between shapes of input actions");
111 : }
112 : }
113 : }
114 : }
115 111 : G::readGridParameters( func.g, action, gridobject, shape );
116 110 : }
117 :
118 : template <class K, class P, class G>
119 105 : void KDEHelper<K,P,G>::setupGridBounds( KDEHelper<K,P,G>& func, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval ) {
120 : // Setup the grid boundaries on first step
121 94 : G::setupGridBounds( func.g, box, gridobject, args, myval );
122 102 : if( gridobject.getGridType()!="fibonacci" ) {
123 : // Set the periodicity of the parameters
124 212 : for(unsigned i=0; i<gridobject.getDimension(); ++i) {
125 121 : P::setArgumentDomain( i, func.kernelsum.params, gridobject.getGridSpacing()[i], gridobject.isPeriodic(i), gridobject.getMin()[i], gridobject.getMax()[i] );
126 : }
127 : }
128 : func.kernelsum.params.update();
129 : // Check if the bandwidth changes during the simulation
130 102 : func.fixed_width = false;
131 102 : std::size_t ndim = gridobject.getDimension();
132 102 : if( K::bandwidthIsConstant( ndim, args ) && K::bandwidthsAllSame( ndim, args ) ) {
133 102 : std::vector<double> myargs( args.size() );
134 477 : for(unsigned j=0; j<args.size(); ++j) {
135 375 : myargs[j] = args[j]->get(0);
136 : }
137 : // Convert the input into kernel parameters, so that the support is computed from the same
138 : // quantities used in transferParamsToKernel when it is recomputed for each kernel
139 : // (NonDiagonalKernelParams stores the inverse of the covariance matrix, not the covariance)
140 102 : std::vector<double> kparams( K::getNumberOfParameters(ndim), 0.0 );
141 : View<double> kpview( kparams.data(), kparams.size() );
142 102 : K::setKernelAndCheckHeight( kpview, ndim, myargs );
143 102 : G::getDiscreteSupport( func.g, func.kernelsum.params, View<const double>( kparams.data() + ndim, K::getNumberOfShapeParameters(ndim)), func.nneigh, gridobject );
144 102 : func.fixed_width = true;
145 : }
146 102 : }
147 :
148 : template <class K, class P, class G>
149 181 : void KDEHelper<K,P,G>::readKernelParameters( std::string& value, ActionWithArguments* action, const std::string& outlab, bool rerequestargs ) {
150 : std::size_t dot = value.find_first_of('.');
151 181 : ActionWithValue* av = action->plumed.getActionSet().selectWithLabel<ActionWithValue*>( value.substr(0,dot) );
152 181 : if( !av ) {
153 133 : std::string matstr, vals = "VALUES=" + value;
154 133 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
155 : std::string nr, nc;
156 14 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
157 14 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
158 42 : matstr = " NROWS=" + nr + " NCOLS=" + nc;
159 : }
160 1446501 : for(unsigned i=1; i<(action->getPntrToArgument(0))->getNumberOfValues(); ++i) {
161 2892736 : vals += "," + value;
162 : }
163 266 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + ": CONSTANT NOLOG " + vals + matstr ), false );
164 266 : value = action->getLabel() + outlab;
165 : } else {
166 : Value* myval;
167 48 : if( dot!=std::string::npos ) {
168 2 : myval = av->copyOutput( value );
169 : } else {
170 46 : if( av->getNumberOfComponents()>1 ) {
171 0 : action->error("problem reading argument " + value );
172 : }
173 46 : myval = av->copyOutput(0);
174 : }
175 48 : if( myval->getRank()==0 ) {
176 : std::string nvals;
177 0 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
178 : std::string nr, nc;
179 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
180 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
181 0 : nvals = nr + "," + nc;
182 : } else {
183 0 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), nvals );
184 : }
185 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + "_ones: ONES SIZE=" + nvals ), false );
186 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + outlab + ": CUSTOM ARG=" + action->getLabel() + outlab + "_ones," + value ), false );
187 0 : value = action->getLabel() + outlab;
188 : }
189 : }
190 181 : if( !rerequestargs ) {
191 : return;
192 : }
193 150 : KDEHelper<K,P,G>::addArgument( value, action );
194 0 : }
195 :
196 : template <class K, class P, class G>
197 242 : void KDEHelper<K,P,G>::addArgument( const std::string& value, ActionWithArguments* action ) {
198 242 : std::size_t dot = value.find_first_of(".");
199 242 : ActionWithValue* av = action->plumed.getActionSet().selectWithLabel<ActionWithValue*>( value.substr(0,dot) );
200 242 : plumed_assert( av );
201 242 : std::vector<Value*> args( action->getArguments() );
202 242 : args.push_back( av->copyOutput(0) );
203 242 : action->requestArguments( args );
204 242 : }
205 :
206 : template <class K, class P, class G>
207 600634 : 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, View<double>& params ) {
208 : // This sets the kernel parameters for the Kth kernel and checks that we want
209 : // to consider it
210 600634 : std::size_t ndim = gridobject.getDimension();
211 600634 : if( !K::setKernelAndCheckHeight( params, ndim, argval ) ) {
212 224647 : return;
213 : }
214 : // If the widths of each kernel are not all the same then get the discrete support
215 375987 : if( updateNeighborsOnEachKernel ) {
216 0 : G::getDiscreteSupport( func.g, func.kernelsum.params, View<const double>(params.data()+ndim,K::getNumberOfShapeParameters(ndim)), func.nneigh, gridobject );
217 : }
218 :
219 : // Now get the grid points for this particular kernel
220 : unsigned num_neigh;
221 375987 : std::vector<unsigned> neighbors;
222 375987 : G::getNeighbors( func.kernelsum.params, View<double>(params.data(),ndim), gridobject, func.nneigh, num_neigh, neighbors );
223 :
224 : // And transfer the neighbor information to the holders
225 57114839 : for(unsigned j=0; j<num_neigh; ++j) {
226 56738853 : plumed_assert( neighbors[j]*nkernels + func.nkernels_per_point_v[neighbors[j]] < func.kernels_for_gridpoint_v.size() );
227 56738853 : func.kernels_for_gridpoint_v[ neighbors[j]*nkernels + func.nkernels_per_point_v[neighbors[j]] ] = kval;
228 56738853 : func.nkernels_per_point_v[ neighbors[j] ]++;
229 : }
230 : }
231 :
232 : template <class K, class P, class G>
233 2600 : void KDEHelper<K,P,G>::transferKernels( KDEHelper<K,P,G>& func, const std::vector<Value*>& args, GridCoordinatesObject& gridobject ) {
234 : // Resize the kernel sum if we need to
235 : // Number of kernels is determined based on sparsity pattern of matrix input as matrix of heights
236 2600 : std::size_t nkernels = args[args.size()-1]->getNumberOfStoredValues();
237 : // And resize the grid counters if we need to
238 2600 : std::size_t ngp = gridobject.getNumberOfPoints();
239 2600 : if( func.nkernels_per_point_v.size()!=ngp ) {
240 102 : func.nkernels_per_point_v.resize( ngp );
241 : }
242 2600 : std::size_t ndim = gridobject.getDimension();
243 : std::size_t nparams_per_kernel = K::getNumberOfParameters(ndim);
244 : std::fill( func.nkernels_per_point_v.begin(), func.nkernels_per_point_v.end(), 0 );
245 2600 : if( func.kernelsum.kernelParams_v.size()!=nkernels*nparams_per_kernel ) {
246 109 : func.kernelsum.kernelParams_v.resize( nkernels*nparams_per_kernel );
247 109 : func.kernels_for_gridpoint_v.resize( ngp*nkernels );
248 : }
249 :
250 2600 : std::vector<double> argval( args.size() );
251 2600 : bool bandwidthsAllSame = false, updateNeighborsOnEachKernel = !func.fixed_width;
252 2600 : if( updateNeighborsOnEachKernel ) {
253 0 : bandwidthsAllSame = K::bandwidthsAllSame( ndim, args );
254 : }
255 2600 : if( args[args.size()-1]->getRank()==2 ) {
256 53 : const unsigned nc = args[args.size()-1]->getShape()[1];
257 53 : const unsigned nrows = args[args.size()-1]->getShape()[0];
258 53 : const unsigned ncs = args[args.size()-1]->getNumberOfColumns();
259 11237 : for(unsigned i=0; i<nrows; ++i) {
260 11184 : unsigned ncols = args[args.size()-1]->getRowLength(i);
261 60646 : for(unsigned k=0; k<args.size(); ++k) {
262 87640 : plumed_massert( args[k]->isConstant() || ncols==args[k]->getRowLength(i), "all input matrices must have same sparsity pattern" );
263 : }
264 448275 : for(unsigned j=0; j<ncols; ++j) {
265 437091 : unsigned jind = args[args.size()-1]->getRowIndex( i, j );
266 1990664 : for(unsigned k=0; k<args.size(); ++k) {
267 1553573 : if( jind==args[k]->getRowIndex( i, j ) ) {
268 1138794 : argval[k] = args[k]->get( i*ncs+ j, false );
269 : } else {
270 414779 : argval[k] = args[k]->get( i*nc + jind );
271 : }
272 : }
273 437091 : View<double> myparams( func.kernelsum.kernelParams_v.data() + nparams_per_kernel*(i*ncs+j), nparams_per_kernel );
274 437091 : KDEHelper<K,P,G>::transferParamsToKernel( argval, func, gridobject, updateNeighborsOnEachKernel, nkernels, i*ncs+j, myparams );
275 437091 : if( bandwidthsAllSame ) {
276 : updateNeighborsOnEachKernel = false;
277 : }
278 : }
279 : }
280 : } else {
281 166089 : for(unsigned i=0; i<nkernels; ++i) {
282 : // Transfer the kernel parameters to local vector of doubles
283 873139 : for(unsigned j=0; j<args.size(); ++j) {
284 709596 : argval[j] = args[j]->get(i,false);
285 : }
286 163543 : View<double> myparams( func.kernelsum.kernelParams_v.data() + nparams_per_kernel*i, nparams_per_kernel );
287 163543 : KDEHelper<K,P,G>::transferParamsToKernel( argval, func, gridobject, updateNeighborsOnEachKernel, nkernels, i, myparams );
288 163542 : if( bandwidthsAllSame ) {
289 : updateNeighborsOnEachKernel = false;
290 : }
291 : }
292 : }
293 : // Get the maximum number of kernels for any given grid point (used for resizing derivatives)
294 2599 : func.maxkernels = 0;
295 586780 : for(unsigned i=0; i<ngp; ++i) {
296 584181 : if( func.nkernels_per_point_v[i]>func.maxkernels ) {
297 3919 : func.maxkernels = func.nkernels_per_point_v[i];
298 : }
299 : }
300 : func.update();
301 2599 : }
302 :
303 : template <class K, class P, class G, typename myPTM=defaultPTM>
304 : class KDE : public ActionWithGrid {
305 : public:
306 : using input_type = KDEHelper<K, P, G>;
307 : using mytype=KDE<K, P, G, myPTM>;
308 : using PTM = typename myPTM::template PTM<mytype>;
309 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
310 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
311 : private:
312 : bool firststep;
313 : /// The parallel task manager
314 : PTM taskmanager;
315 : GridCoordinatesObject gridobject;
316 : public:
317 : static void registerKeywords( Keywords& keys );
318 : explicit KDE(const ActionOptions&ao);
319 : std::vector<std::string> getGridCoordinateNames() const override ;
320 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
321 : unsigned getNumberOfDerivatives() override;
322 : int checkTaskIsActive( const unsigned& itask ) const override ;
323 : void turnOnDerivatives() override ;
324 : void prepare() override ;
325 : void calculate() override ;
326 : void getArgumentInputData( std::vector<double>& inputdata ) const override ;
327 : void getArgumentInputData( std::vector<float>& inputdata ) const override ;
328 : static void performTask( std::size_t task_index,
329 : const KDEHelper<K, P, G>& actiondata,
330 : ParallelActionsInput& input,
331 : ParallelActionsOutput& output );
332 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
333 : static int getNumberOfValuesPerTask( std::size_t task_index,
334 : const KDEHelper<K, P, G>& actiondata );
335 : static void getForceIndices( std::size_t task_index,
336 : std::size_t colno,
337 : std::size_t ntotal_force,
338 : const KDEHelper<K, P, G>& actiondata,
339 : const ParallelActionsInput& input,
340 : ForceIndexHolder force_indices );
341 : };
342 :
343 : template <class K, class P, class G, typename myPTM>
344 443 : void KDE<K,P,G,myPTM>::registerKeywords( Keywords& keys ) {
345 443 : ActionWithGrid::registerKeywords( keys );
346 886 : keys.addInputKeyword("compulsory","ARG","scalar/vector/matrix","the label for the value that should be used to construct the histogram");
347 443 : KDEHelper<K,P,G>::registerKeywords( keys );
348 : // Keywords for spherical KDE
349 443 : keys.add("hidden","MASKED_INPUT_ALLOWED","turns on that you are allowed to use masked inputs ");
350 886 : keys.setValueDescription("grid","a function on a grid that was obtained by doing a Kernel Density Estimation using the input arguments");
351 443 : if( keys.getDisplayName()!="SPHERICAL_KDE" ) {
352 852 : keys.setDisplayName("KDE");
353 : }
354 443 : PTM::registerKeywords( keys );
355 443 : }
356 :
357 : template <class K, class P, class G, typename myPTM>
358 112 : KDE<K,P,G,myPTM>::KDE(const ActionOptions&ao):
359 : Action(ao),
360 : ActionWithGrid(ao),
361 112 : firststep(true),
362 112 : taskmanager(this) {
363 :
364 112 : std::vector<std::size_t> shape( getNumberOfArguments() );
365 112 : unsigned numberOfKernels=getPntrToArgument(0)->getNumberOfValues();
366 172 : for(unsigned i=1; i<shape.size(); ++i) {
367 60 : if( numberOfKernels!=getPntrToArgument(i)->getNumberOfValues() ) {
368 0 : error("mismatch between numbers of values in input arguments");
369 : }
370 : }
371 :
372 112 : function::FunctionOptions foptions;
373 112 : KDEHelper<K,P,G>::read( taskmanager.getActionInput(), this, getArguments(), gridobject, shape, foptions );
374 110 : addValueWithDerivatives( shape );
375 110 : setNotPeriodic();
376 110 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
377 120 : }
378 :
379 : template <class K, class P, class G, typename myPTM>
380 8063 : unsigned KDE<K,P,G,myPTM>::getNumberOfDerivatives() {
381 8063 : return gridobject.getDimension();
382 : }
383 :
384 : template <class K, class P, class G, typename myPTM>
385 155 : std::vector<std::string> KDE<K,P,G,myPTM>::getGridCoordinateNames() const {
386 155 : std::vector<std::string> names( gridobject.getDimension() );
387 386 : for(unsigned i=0; i<names.size(); ++i) {
388 : names[i] = getPntrToArgument(i)->getName();
389 : }
390 155 : return names;
391 0 : }
392 :
393 : template <class K, class P, class G, typename myPTM>
394 280185 : const GridCoordinatesObject& KDE<K,P,G,myPTM>::getGridCoordinatesObject() const {
395 280185 : return gridobject;
396 : }
397 :
398 : template <class K, class P, class G, typename myPTM>
399 584181 : int KDE<K,P,G,myPTM>::checkTaskIsActive( const unsigned& itask ) const {
400 584181 : if( taskmanager.getActionInput().nkernels_per_point_v[itask]>0 ) {
401 183368 : return 1;
402 : }
403 : return -1;
404 : }
405 :
406 : template <class K, class P, class G, typename myPTM>
407 64 : void KDE<K,P,G,myPTM>::turnOnDerivatives() {
408 64 : if( getName().find("DISCRETE")!=std::string::npos ) {
409 0 : error("KDE with discrete kernels is not differentiable");
410 : }
411 64 : ActionWithValue::turnOnDerivatives();
412 64 : }
413 :
414 : template <class K, class P, class G, typename myPTM>
415 2603 : void KDE<K,P,G,myPTM>::prepare() {
416 2603 : ActionWithVector::prepare();
417 2603 : std::size_t nkernels = getPntrToArgument(0)->getNumberOfValues();
418 8101 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
419 : Value* myarg = getPntrToArgument(i);
420 5498 : if( myarg->getNumberOfValues()!=nkernels ) {
421 1 : if( myarg->isConstant() && myarg->getNumberOfValues()==1 ) {
422 1 : myarg->reshapeConstantValue( getPntrToArgument(0)->getShape() );
423 : } else {
424 0 : plumed_merror("found mismatched numbers of arguments in input");
425 : }
426 : }
427 : }
428 2603 : }
429 :
430 : template <class K, class P, class G, typename myPTM>
431 2603 : void KDE<K,P,G,myPTM>::calculate() {
432 2603 : if( firststep ) {
433 105 : Tensor box;
434 105 : PbcAction* bv = plumed.getActionSet().template selectWithLabel<PbcAction*>("Box");
435 105 : if( bv ) {
436 104 : box = bv->getPbc().getBox();
437 : }
438 105 : KDEHelper<K,P,G>::setupGridBounds( taskmanager.getActionInput(), box, gridobject, getArguments(), getPntrToComponent(0) );
439 102 : firststep=false;
440 : }
441 2600 : KDEHelper<K,P,G>::transferKernels( taskmanager.getActionInput(), getArguments(), gridobject );
442 2599 : taskmanager.setupParallelTaskManager( getNumberOfArguments()*taskmanager.getActionInput().maxkernels, getNumberOfForceDerivatives() );
443 2599 : taskmanager.runAllTasks();
444 2599 : }
445 :
446 : template <class K, class P, class G, typename myPTM>
447 2599 : void KDE<K,P,G,myPTM>::getArgumentInputData( std::vector<double>& inputdata ) const {
448 2599 : std::size_t ndim = gridobject.getDimension();
449 2599 : std::size_t nstored = getConstPntrToComponent(0)->getNumberOfStoredValues();
450 2599 : std::vector<double> pos( ndim );
451 2599 : if( inputdata.size()!=nstored*ndim ) {
452 101 : inputdata.resize( ndim*nstored );
453 : }
454 :
455 586780 : for(unsigned i=0; i<nstored; ++i) {
456 584181 : gridobject.getGridPointCoordinates( i, pos );
457 1939989 : for(unsigned j=0; j<ndim; ++j) {
458 1355808 : inputdata[ i*ndim + j ] = pos[j];
459 : }
460 : }
461 2599 : }
462 :
463 : template <class K, class P, class G, typename myPTM>
464 0 : void KDE<K,P,G,myPTM>::getArgumentInputData( std::vector<float>& inputdata ) const {
465 0 : std::size_t ndim = gridobject.getDimension();
466 0 : std::size_t nstored = getConstPntrToComponent(0)->getNumberOfStoredValues();
467 0 : std::vector<double> pos( ndim );
468 0 : if( inputdata.size()!=nstored*ndim ) {
469 0 : inputdata.resize( ndim*nstored );
470 : }
471 :
472 0 : for(unsigned i=0; i<nstored; ++i) {
473 0 : gridobject.getGridPointCoordinates( i, pos );
474 0 : for(unsigned j=0; j<ndim; ++j) {
475 0 : inputdata[ i*ndim + j ] = pos[j];
476 : }
477 : }
478 0 : }
479 :
480 : template <class K, class P, class G, typename myPTM>
481 190619 : void KDE<K,P,G,myPTM>::performTask( std::size_t task_index,
482 : const KDEHelper<K, P, G>& actiondata,
483 : ParallelActionsInput& input,
484 : ParallelActionsOutput& output ) {
485 190619 : std::size_t ndim = actiondata.ndim;
486 381238 : SumOfKernels<K,P>::calc( View<const std::size_t>( actiondata.kernels_for_gridpoint + task_index*input.argstarts[1], actiondata.nkernels_per_point[task_index] ),
487 190619 : actiondata.kernelsum,
488 190619 : View<const double>(input.inputdata + task_index*ndim,ndim),
489 : View<double>(output.values.data(), 1),
490 : View<double>(output.values.data()+1, ndim),
491 190619 : View<double>(output.derivatives.data(),actiondata.maxkernels*input.nargs) );
492 :
493 190619 : }
494 :
495 : template <class K, class P, class G, typename myPTM>
496 617 : void KDE<K,P,G,myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
497 617 : taskmanager.applyForces( outforces );
498 617 : }
499 :
500 : template <class K, class P, class G, typename myPTM>
501 7329 : int KDE<K,P,G,myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
502 : const KDEHelper<K, P, G>& actiondata ) {
503 7329 : return 1;
504 : }
505 :
506 : template <class K, class P, class G, typename myPTM>
507 7239 : void KDE<K,P,G,myPTM>::getForceIndices( std::size_t task_index,
508 : std::size_t colno,
509 : std::size_t ntotal_force,
510 : const KDEHelper<K, P, G>& actiondata,
511 : const ParallelActionsInput& input,
512 : ForceIndexHolder force_indices ) {
513 7239 : force_indices.threadsafe_derivatives_end[0] = 0;
514 : std::size_t nparams = K::getNumberOfParameters( actiondata.ndim );
515 7239 : View<const std::size_t> kernellist( actiondata.kernels_for_gridpoint + task_index*input.argstarts[1], actiondata.nkernels_per_point[task_index] );
516 1051320 : for(unsigned i=0; i<kernellist.size(); ++i) {
517 5555892 : for(unsigned j=0; j<nparams; ++j) {
518 4511811 : force_indices.indices[0][i*nparams+j] = input.argstarts[j] + kernellist[i];
519 : }
520 : }
521 7239 : force_indices.tot_indices[0] = kernellist.size()*nparams;
522 7239 : }
523 :
524 : }
525 : }
526 : #endif
|
