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 : #include "ActionWithGrid.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PbcAction.h"
27 : #include "tools/HistogramBead.h"
28 : #include "tools/SwitchingFunction.h"
29 : #include "tools/Matrix.h"
30 :
31 : //+PLUMEDOC ANALYSIS KDE
32 : /*
33 : Create a histogram from the input scalar/vector/matrix using KDE
34 :
35 : \par Examples
36 :
37 :
38 : */
39 : //+ENDPLUMEDOC
40 :
41 : //+PLUMEDOC ANALYSIS SPHERICAL_KDE
42 : /*
43 : Create a histogram from the input scalar/vector/matrix using SPHERICAL_KDE
44 :
45 : \par Examples
46 :
47 :
48 : */
49 : //+ENDPLUMEDOC
50 :
51 : namespace PLMD {
52 : namespace gridtools {
53 :
54 : class KDE : public ActionWithGrid {
55 : private:
56 : double hh;
57 : bool hasheight;
58 : bool ignore_out_of_bounds, fixed_width;
59 : double dp2cutoff;
60 : std::string kerneltype;
61 : GridCoordinatesObject gridobject;
62 : std::vector<std::string> gmin, gmax;
63 : std::vector<double> center;
64 : std::vector<double> gspacing;
65 : unsigned num_neigh, bwargno;
66 : std::vector<Value> grid_diff_value;
67 : std::vector<unsigned> nbin, nneigh, neighbors;
68 : unsigned numberOfKernels, nbins;
69 : SwitchingFunction switchingFunction;
70 : double von_misses_concentration, von_misses_norm;
71 : void setupNeighborsVector();
72 : void retrieveArgumentsAndHeight( const MultiValue& myvals, std::vector<double>& args, double& height ) const ;
73 : double evaluateKernel( const std::vector<double>& gpoint, const std::vector<double>& args, const double& height, std::vector<double>& der ) const ;
74 : void setupHistogramBeads( std::vector<HistogramBead>& bead ) const ;
75 : double evaluateBeadValue( std::vector<HistogramBead>& bead, const std::vector<double>& gpoint, const std::vector<double>& args, const double& height, std::vector<double>& der ) const ;
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit KDE(const ActionOptions&ao);
79 : std::vector<std::string> getGridCoordinateNames() const override ;
80 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
81 : unsigned getNumberOfDerivatives() override;
82 : void setupOnFirstStep( const bool incalc ) override ;
83 : void getNumberOfTasks( unsigned& ntasks ) override ;
84 : void areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) override ;
85 : int checkTaskStatus( const unsigned& taskno, int& flag ) const override ;
86 : void performTask( const unsigned& current, MultiValue& myvals ) const override ;
87 : void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
88 : const unsigned& bufstart, std::vector<double>& buffer ) const override ;
89 : void updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const override ;
90 : void gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const override ;
91 : };
92 :
93 : PLUMED_REGISTER_ACTION(KDE,"KDE")
94 : PLUMED_REGISTER_ACTION(KDE,"SPHERICAL_KDE")
95 :
96 274 : void KDE::registerKeywords( Keywords& keys ) {
97 274 : ActionWithGrid::registerKeywords( keys ); keys.use("ARG");
98 548 : keys.add("optional","HEIGHTS","this keyword takes the label of an action that calculates a vector of values. The elements of this vector "
99 : "are used as weights for the Gaussians.");
100 548 : keys.add("optional","VOLUMES","this keyword take the label of an action that calculates a vector of values. The elements of this vector "
101 : "divided by the volume of the Gaussian are used as weights for the Gaussians");
102 : // Keywords for KDE
103 548 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
104 548 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
105 548 : keys.add("optional","BANDWIDTH","the bandwidths for kernel density esimtation");
106 548 : keys.add("compulsory","METRIC","the inverse covariance to use for the kernels that are added to the grid");
107 548 : keys.add("compulsory","CUTOFF","6.25","the cutoff at which to stop evaluating the kernel functions is set equal to sqrt(2*x)*bandwidth in each direction where x is this number");
108 548 : keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using. More details on the kernels available "
109 : "in plumed plumed can be found in \\ref kernelfunctions.");
110 548 : keys.add("optional","GRID_BIN","the number of bins for the grid");
111 548 : keys.addFlag("IGNORE_IF_OUT_OF_RANGE",false,"if a kernel is outside of the range of the grid it is safe to ignore");
112 548 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
113 : // Keywords for spherical KDE
114 548 : keys.add("compulsory","CONCENTRATION","the concentration parameter for Von Mises-Fisher distributions (only required for SPHERICAL_KDE)");
115 274 : keys.setValueDescription("a function on a grid that was obtained by doing a Kernel Density Estimation using the input arguments");
116 274 : }
117 :
118 149 : KDE::KDE(const ActionOptions&ao):
119 : Action(ao),
120 : ActionWithGrid(ao),
121 149 : hasheight(false),
122 149 : fixed_width(false)
123 : {
124 149 : std::vector<unsigned> shape( getNumberOfArguments() ); center.resize( getNumberOfArguments() );
125 149 : numberOfKernels=getPntrToArgument(0)->getNumberOfValues();
126 199 : for(unsigned i=1; i<shape.size(); ++i) {
127 50 : if( numberOfKernels!=getPntrToArgument(i)->getNumberOfValues() ) error("mismatch between numbers of values in input arguments");
128 : }
129 :
130 298 : bool weights_are_volumes=true; std::vector<std::string> weight_str; parseVector("VOLUMES",weight_str);
131 149 : if( weight_str.size()==0 ) {
132 146 : parseVector("HEIGHTS",weight_str);
133 73 : if( weight_str.size()>0 ) weights_are_volumes=false;
134 : }
135 149 : hasheight=(weight_str.size()==1);
136 149 : if( weight_str.size()>1 ) error("only one scalar/vector/matrix should be input to HEIGHTS");
137 :
138 149 : if( getName()=="KDE" ) {
139 278 : parse("KERNEL",kerneltype);
140 139 : if( kerneltype!="DISCRETE" ) {
141 254 : std::string bandwidth; std::vector<std::string> bwidths; parseVector("BANDWIDTH",bwidths);
142 127 : if( bwidths.size()>0 ) {
143 127 : std::string band="VALUES=" + bwidths[0];
144 284 : for(unsigned i=0; i<bwidths.size(); ++i) {
145 187 : if( i>0 ) band += "," + bwidths[i];
146 : }
147 254 : plumed.readInputLine( getLabel() + "_sigma: CONSTANT " + band );
148 254 : plumed.readInputLine( getLabel() + "_cov: CUSTOM ARG=" + getLabel() + "_sigma FUNC=x*x PERIODIC=NO" );
149 254 : plumed.readInputLine( getLabel() + "_icov: CUSTOM ARG=" + getLabel() + "_cov FUNC=1/x PERIODIC=NO" );
150 254 : bandwidth = getLabel() + "_icov";
151 :
152 254 : if( (kerneltype=="gaussian" || kerneltype=="GAUSSIAN") && weights_are_volumes ) {
153 105 : std::string pstr; Tools::convert( sqrt(pow(2*pi,bwidths.size())), pstr );
154 210 : plumed.readInputLine( getLabel() + "_bwprod: PRODUCT ARG=" + getLabel() + "_cov");
155 210 : plumed.readInputLine( getLabel() + "_vol: CUSTOM ARG=" + getLabel() + "_bwprod FUNC=(sqrt(x)*" + pstr + ") PERIODIC=NO");
156 181 : if( hasheight ) plumed.readInputLine( getLabel() + "_height: CUSTOM ARG=" + weight_str[0] + "," + getLabel() + "_vol FUNC=x/y PERIODIC=NO");
157 58 : else plumed.readInputLine( getLabel() + "_height: CUSTOM ARG=" + getLabel() + "_vol FUNC=1/x PERIODIC=NO");
158 210 : hasheight=true; weight_str.resize(1); weight_str[0] = getLabel() + "_height";
159 : }
160 0 : } else parse("METRIC",bandwidth);
161 127 : weight_str.push_back( bandwidth );
162 127 : }
163 : }
164 149 : if( weight_str.size()>0 ) {
165 145 : std::vector<Value*> weight_args; ActionWithArguments::interpretArgumentList( weight_str, plumed.getActionSet(), this, weight_args );
166 145 : std::vector<Value*> args( getArguments() ); args.push_back( weight_args[0] );
167 145 : if( hasheight && weight_args[0]->getNumberOfValues()>1 && numberOfKernels!=weight_args[0]->getNumberOfValues() ) error("mismatch between numbers of values in input arguments and HEIGHTS");
168 :
169 145 : if( weight_str.size()==2 ) {
170 112 : log.printf(" quantities used for weights are : %s \n", weight_str[0].c_str() );
171 112 : args.push_back( weight_args[1] );
172 112 : if( weight_args[1]->getRank()==1 && weight_args[1]->getNumberOfValues()!=shape.size() ) error("size of bandwidth vector is incorrect");
173 112 : if( weight_args[1]->getRank()>2 ) error("bandwidths cannot have rank greater than 2");
174 112 : bwargno=args.size()-1; log.printf(" bandwidths are taken from : %s \n", weight_str[1].c_str() );
175 33 : } else if( !hasheight ) {
176 15 : if( weight_args[0]->getRank()==1 && weight_args[0]->getNumberOfValues()!=shape.size() ) error("size of bandwidth vector is incorrect");
177 15 : if( weight_args[0]->getRank()>2 ) error("bandwidths cannot have rank greater than 2");
178 15 : bwargno=args.size()-1; log.printf(" bandwidths are taken from : %s \n", weight_str[0].c_str() );
179 18 : } else if ( weight_str.size()==1 ) {
180 18 : log.printf(" quantities used for weights are : %s \n", weight_str[0].c_str() );
181 0 : } else error("only one scalar/vector/matrix should be input to HEIGHTS");
182 145 : requestArguments( args );
183 : }
184 :
185 149 : if( getName()=="KDE" ) {
186 139 : bool hasauto=false; gmin.resize( shape.size() ); gmax.resize( shape.size() );
187 278 : parseVector("GRID_MIN",gmin); parseVector("GRID_MAX",gmax);
188 308 : for(unsigned i=0; i<gmin.size(); ++i) {
189 169 : if( gmin[i]=="auto" ) {
190 52 : log.printf(" for %dth coordinate min and max are set from cell directions \n", (i+1) );
191 : hasauto=true; // We need to do a preparation step to set the grid from the box size
192 52 : if( gmax[i]!="auto" ) error("if gmin is set from box vectors gmax must also be set in the same way");
193 52 : if( getPntrToArgument(i)->isPeriodic() ) {
194 0 : if( gmin[i]=="auto" ) getPntrToArgument(i)->getDomain( gmin[i], gmax[i] );
195 : else {
196 0 : std::string str_min, str_max; getPntrToArgument(i)->getDomain( str_min, str_max );
197 0 : if( str_min!=gmin[i] || str_max!=gmax[i] ) error("all periodic arguments should have the same domain");
198 : }
199 52 : } else if( getPntrToArgument(i)->getName().find(".")!=std::string::npos ) {
200 52 : std::size_t dot = getPntrToArgument(i)->getName().find_first_of(".");
201 52 : std::string name = getPntrToArgument(i)->getName().substr(dot+1);
202 90 : if( name!="x" && name!="y" && name!="z" ) {
203 0 : error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
204 : }
205 : } else {
206 0 : error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
207 : }
208 : } else {
209 117 : log.printf(" for %dth coordinate min is set to %s and max is set to %s \n", (i+1), gmin[i].c_str(), gmax[i].c_str() );
210 : }
211 : }
212 139 : if( hasauto && gmin.size()>3 ) error("can only set GRID_MIN and GRID_MAX automatically if components of distance are used in input");
213 :
214 417 : parseVector("GRID_BIN",nbin); parseVector("GRID_SPACING",gspacing); parse("CUTOFF",dp2cutoff);
215 260 : if( kerneltype.find("bin")==std::string::npos && kerneltype!="DISCRETE" ) {
216 218 : std::string errors; switchingFunction.set( kerneltype + " R_0=1.0 NOSTRETCH", errors );
217 109 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
218 : }
219 :
220 139 : if( nbin.size()!=shape.size() && gspacing.size()!=shape.size() ) error("GRID_BIN or GRID_SPACING must be set");
221 : // Create a value
222 139 : std::vector<bool> ipbc( shape.size() );
223 308 : for(unsigned i=0; i<shape.size(); ++i) {
224 324 : if( getPntrToArgument( i )->isPeriodic() || gmin[i]=="auto" ) ipbc[i]=true;
225 : else ipbc[i]=false;
226 : }
227 278 : gridobject.setup( "flat", ipbc, 0, 0.0 );
228 : } else {
229 10 : if( shape.size()!=3 ) error("should have three coordinates in input to this action");
230 :
231 20 : parse("GRID_BIN",nbins); log.printf(" setting number of bins to %d \n", nbins );
232 10 : parse("CONCENTRATION",von_misses_concentration); fixed_width=true;
233 10 : von_misses_norm = von_misses_concentration / ( 4*pi*sinh( von_misses_concentration ) );
234 10 : log.printf(" setting concentration parameter to %f \n", von_misses_concentration );
235 :
236 : // Create a value
237 10 : std::vector<bool> ipbc( shape.size(), false );
238 10 : double fib_cutoff = std::log( epsilon / von_misses_norm ) / von_misses_concentration;
239 20 : gridobject.setup( "fibonacci", ipbc, nbins, fib_cutoff ); checkRead();
240 :
241 : // Setup the grid
242 10 : shape[0]=nbins; shape[1]=shape[2]=1;
243 : }
244 149 : parseFlag("IGNORE_IF_OUT_OF_RANGE",ignore_out_of_bounds);
245 149 : if( ignore_out_of_bounds ) log.printf(" ignoring kernels that are outside of grid \n");
246 149 : addValueWithDerivatives( shape ); setNotPeriodic();
247 149 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
248 : // Make sure we store all the arguments
249 605 : for(unsigned i=0; i<getNumberOfArguments(); ++i) getPntrToArgument(i)->buildDataStore();
250 : // Check for task reduction
251 149 : updateTaskListReductionStatus(); setupOnFirstStep( false );
252 298 : }
253 :
254 295 : void KDE::setupOnFirstStep( const bool incalc ) {
255 295 : if( getName()=="SPHERICAL_KDE" ) return ;
256 :
257 606 : for(unsigned i=0; i<getNumberOfDerivatives(); ++i) {
258 331 : if( gmin[i]=="auto" && incalc ) {
259 : double lcoord, ucoord;
260 92 : PbcAction* bv = plumed.getActionSet().selectWithLabel<PbcAction*>("Box"); Tensor box( bv->getPbc().getBox() );
261 46 : std::size_t dot = getPntrToArgument(i)->getName().find_first_of(".");
262 46 : std::string name = getPntrToArgument(i)->getName().substr(dot+1);
263 46 : if( name=="x" ) { lcoord=-0.5*box(0,0); ucoord=0.5*box(0,0); }
264 22 : else if( name=="y" ) { lcoord=-0.5*box(1,1); ucoord=0.5*box(1,1); }
265 10 : else if( name=="z" ) { lcoord=-0.5*box(2,2); ucoord=0.5*box(2,2); }
266 0 : else plumed_error();
267 : // And convert to strings for bin and bmax
268 46 : Tools::convert( lcoord, gmin[i] ); Tools::convert( ucoord, gmax[i] );
269 : }
270 331 : if( incalc ) {
271 324 : grid_diff_value.push_back( Value() );
272 162 : if( gridobject.isPeriodic(i) ) grid_diff_value[i].setDomain( gmin[i], gmax[i] );
273 103 : else grid_diff_value[i].setNotPeriodic();
274 : }
275 : }
276 : // And setup the grid object
277 275 : gridobject.setBounds( gmin, gmax, nbin, gspacing );
278 275 : std::vector<unsigned> shape( gridobject.getNbin(true) );
279 275 : getPntrToComponent(0)->setShape( shape );
280 867 : bool hasauto=false; for(unsigned i=0; i<gmin.size(); ++i) { if(gmin[i]=="auto" || gmax[i]=="auto" ) { hasauto=true; break; } }
281 : // And setup the neighbors
282 275 : if( !hasauto && kerneltype!="DISCRETE" && getPntrToArgument(bwargno)->isConstant() ) {
283 214 : fixed_width=true; setupNeighborsVector();
284 : }
285 : }
286 :
287 237 : void KDE::setupNeighborsVector() {
288 237 : if( kerneltype!="DISCRETE" ) {
289 214 : std::vector<double> support(gmin.size(),0); nneigh.resize( gmin.size() );
290 214 : if( kerneltype.find("bin")!=std::string::npos ) {
291 18 : std::size_t dd = kerneltype.find("-bin");
292 18 : HistogramBead bead; bead.setKernelType( kerneltype.substr(0,dd) );
293 18 : Value* bw_arg=getPntrToArgument(bwargno);
294 18 : if( bw_arg->getRank()<2 ) {
295 36 : for(unsigned i=0; i<support.size(); ++i) {
296 18 : bead.set( 0, gridobject.getGridSpacing()[i], 1./sqrt(bw_arg->get(i)) );
297 18 : support[i] = bead.getCutoff(); nneigh[i] = static_cast<unsigned>( ceil( support[i]/gridobject.getGridSpacing()[i] ));
298 : }
299 0 : } else plumed_error();
300 : } else {
301 196 : Value* bw_arg=getPntrToArgument(bwargno);
302 196 : if( bw_arg->getRank()<2 ) {
303 432 : for(unsigned i=0; i<support.size(); ++i) {
304 236 : support[i] = sqrt(2.0*dp2cutoff)*(1.0/sqrt(bw_arg->get(i)));
305 236 : nneigh[i] = static_cast<unsigned>( ceil( support[i] / gridobject.getGridSpacing()[i] ) );
306 : }
307 0 : } else if( bw_arg->getRank()==2 ) {
308 0 : Matrix<double> metric(support.size(),support.size()); unsigned k=0;
309 0 : for(unsigned i=0; i<support.size(); ++i) {
310 0 : for(unsigned j=0; j<support.size(); ++j) { metric(i,j)=bw_arg->get(k); k++; }
311 : }
312 0 : Matrix<double> myautovec(support.size(),support.size()); std::vector<double> myautoval(support.size());
313 0 : diagMat(metric,myautoval,myautovec); double maxautoval=1/myautoval[0]; unsigned ind_maxautoval=0;
314 0 : for(unsigned i=1; i<support.size(); i++) {
315 0 : double neweig=1/myautoval[i]; if(neweig>maxautoval) { maxautoval=neweig; ind_maxautoval=i; }
316 : }
317 0 : for(unsigned i=0; i<support.size(); i++) {
318 0 : support[i] = sqrt(2.0*dp2cutoff)*fabs(sqrt(maxautoval)*myautovec(i,ind_maxautoval));
319 0 : nneigh[i] = static_cast<unsigned>( ceil( support[i] / gridobject.getGridSpacing()[i] ) );
320 : }
321 0 : } else plumed_error();
322 : }
323 468 : for(unsigned i=0; i<gridobject.getDimension(); ++i) {
324 254 : double fmax, fmin; Tools::convert( gridobject.getMin()[i], fmin ); Tools::convert( gridobject.getMax()[i], fmax );
325 254 : if( gridobject.isPeriodic(i) && 2*support[i]>(fmax-fmin) ) error("bandwidth is too large for periodic grid");
326 : }
327 : }
328 237 : }
329 :
330 1065 : unsigned KDE::getNumberOfDerivatives() {
331 1065 : return gridobject.getDimension();
332 : }
333 :
334 125 : std::vector<std::string> KDE::getGridCoordinateNames() const {
335 125 : std::vector<std::string> names( gridobject.getDimension() );
336 301 : for(unsigned i=0; i<names.size(); ++i) names[i] = getPntrToArgument(i)->getName();
337 125 : return names;
338 0 : }
339 :
340 6222 : const GridCoordinatesObject& KDE::getGridCoordinatesObject() const {
341 6222 : return gridobject;
342 : }
343 :
344 149 : void KDE::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
345 149 : if( numberOfKernels==1 || (hasheight && getPntrToArgument(gridobject.getDimension())->getRank()>0) ) task_reducing_actions.push_back(this);
346 149 : }
347 :
348 2651 : void KDE::getNumberOfTasks( unsigned& ntasks ) {
349 2651 : if( !fixed_width ) setupNeighborsVector();
350 2651 : ntasks = numberOfKernels = getPntrToArgument(0)->getNumberOfValues();
351 2651 : if( numberOfKernels>1 ) return;
352 :
353 132 : hh = 1.0; if( hasheight ) hh = getPntrToArgument(gridobject.getDimension())->get();
354 309 : for(unsigned i=0; i<center.size(); ++i) center[i]=getPntrToArgument(i)->get();
355 132 : if( !ignore_out_of_bounds && !gridobject.inbounds( center ) ) {
356 : //if( fabs(height)>epsilon ) warning("bounds are possibly set too small as hills with substantial heights are being ignored");
357 : return;
358 : }
359 132 : if( kerneltype=="DISCRETE" ) {
360 12 : num_neigh=1; neighbors.resize(1);
361 24 : for(unsigned i=0; i<center.size(); ++i) center[i] += 0.5*gridobject.getGridSpacing()[i];
362 12 : neighbors[0]=gridobject.getIndex( center );
363 120 : } else gridobject.getNeighbors( center, nneigh, num_neigh, neighbors );
364 132 : ntasks = getPntrToComponent(0)->getNumberOfValues();
365 132 : return;
366 : }
367 :
368 1500435 : int KDE::checkTaskStatus( const unsigned& taskno, int& flag ) const {
369 1500435 : if( numberOfKernels>1 ) {
370 1400374 : if( hasheight && getPntrToArgument(gridobject.getDimension())->getRank()>0
371 2823898 : && fabs(getPntrToArgument(gridobject.getDimension())->get(taskno))<epsilon ) return 0;
372 185408 : return 1;
373 : }
374 19394629 : for(unsigned i=0; i<num_neigh; ++i) {
375 19339045 : if( taskno==neighbors[i] ) return 1;
376 : }
377 : return 0;
378 : }
379 :
380 309567 : void KDE::performTask( const unsigned& current, MultiValue& myvals ) const {
381 309567 : if( numberOfKernels==1 ) {
382 21277 : double newval; std::vector<double> args( gridobject.getDimension() ), der( gridobject.getDimension() );
383 21277 : unsigned valout = getConstPntrToComponent(0)->getPositionInStream();
384 21277 : gridobject.getGridPointCoordinates( current, args );
385 21277 : if( getName()=="KDE" ) {
386 21277 : if( kerneltype=="DISCRETE" ) {
387 : newval = 1.0;
388 21265 : } else if( kerneltype.find("bin")!=std::string::npos ) {
389 0 : double val=hh; std::size_t dd = kerneltype.find("-bin");
390 0 : HistogramBead bead; bead.setKernelType( kerneltype.substr(0,dd) ); Value* bw_arg=getPntrToArgument(bwargno);
391 0 : for(unsigned j=0; j<args.size(); ++j) {
392 0 : if( gridobject.isPeriodic(j) ) {
393 0 : double lcoord, ucoord; Tools::convert( gmin[j], lcoord );
394 0 : Tools::convert( gmax[j], ucoord ); bead.isPeriodic( lcoord, ucoord );
395 : } else bead.isNotPeriodic();
396 0 : if( bw_arg->getRank()<2 ) bead.set( args[j], args[j]+gridobject.getGridSpacing()[j], 1/sqrt(bw_arg->get(j)) );
397 0 : else if( bw_arg->getRank()==2 ) plumed_error();
398 0 : double contr = bead.calculateWithCutoff( args[j], der[j] );
399 0 : val = val*contr; der[j] = der[j] / contr;
400 : }
401 0 : for(unsigned j=0; j<args.size(); ++j) der[j] *= val; newval=val;
402 : } else {
403 21265 : newval = evaluateKernel( args, center, hh, der );
404 : }
405 : } else {
406 0 : double dot=0; for(unsigned i=0; i<der.size(); ++i) { dot += args[i]*center[i]; }
407 0 : newval = hh*von_misses_norm*exp( von_misses_concentration*dot );
408 0 : for(unsigned i=0; i<der.size(); ++i) der[i] = von_misses_concentration*newval*args[i];
409 : }
410 21277 : myvals.setValue( valout, newval );
411 78047 : for(unsigned i=0; i<der.size(); ++i) { myvals.addDerivative( valout, i, der[i] ); myvals.updateIndex( valout, i ); }
412 : }
413 309567 : }
414 :
415 288290 : void KDE::retrieveArgumentsAndHeight( const MultiValue& myvals, std::vector<double>& args, double& height ) const {
416 712540 : height=1.0; for(unsigned i=0; i<args.size(); ++i) args[i]=getPntrToArgument(i)->get( myvals.getTaskIndex() );
417 288290 : if( hasheight && getPntrToArgument(args.size())->getRank()==0 ) height = getPntrToArgument( args.size() )->get();
418 272358 : else if( hasheight ) height = getPntrToArgument( args.size() )->get( myvals.getTaskIndex() );
419 288290 : }
420 :
421 38119671 : double KDE::evaluateKernel( const std::vector<double>& gpoint, const std::vector<double>& args, const double& height, std::vector<double>& der ) const {
422 38119671 : double r2=0, hval = height; Value* bw_arg=getPntrToArgument(bwargno);
423 38119671 : if( bw_arg->getRank()<2 ) {
424 151428182 : for(unsigned j=0; j<der.size(); ++j) {
425 113308511 : double tmp = -grid_diff_value[j].difference( gpoint[j], args[j] );
426 113308511 : der[j] = tmp*bw_arg->get(j); r2 += tmp*der[j];
427 : }
428 0 : } else if( bw_arg->getRank()==2 ) {
429 0 : for(unsigned j=0; j<der.size(); ++j) {
430 0 : der[j]=0; double dp_j, dp_k; dp_j = -grid_diff_value[j].difference( gpoint[j], args[j] );
431 0 : for(unsigned k=0; k<der.size(); ++k ) {
432 0 : if(j==k) dp_k = dp_j;
433 0 : else dp_k = -grid_diff_value[k].difference( gpoint[k], args[k] );
434 0 : der[j] += bw_arg->get(j*der.size()+k)*dp_k; r2 += dp_j*dp_k*bw_arg->get(j*der.size()+k);
435 : }
436 : }
437 0 : } else plumed_error();
438 38119671 : double dval, val=hval*switchingFunction.calculateSqr( r2, dval );
439 151428182 : dval *= hval; for(unsigned j=0; j<der.size(); ++j) der[j] *= dval;
440 38119671 : return val;
441 : }
442 :
443 133400 : void KDE::setupHistogramBeads( std::vector<HistogramBead>& bead ) const {
444 133400 : std::size_t dd = kerneltype.find("-bin");
445 133400 : std::string ktype=kerneltype.substr(0,dd);
446 266800 : for(unsigned j=0; j<bead.size(); ++j) {
447 133400 : bead[j].setKernelType( ktype );
448 133400 : if( gridobject.isPeriodic(j) ) {
449 133400 : double lcoord, ucoord; Tools::convert( gmin[j], lcoord );
450 133400 : Tools::convert( gmax[j], ucoord ); bead[j].isPeriodic( lcoord, ucoord );
451 : } else bead[j].isNotPeriodic();
452 : }
453 133400 : }
454 :
455 533600 : double KDE::evaluateBeadValue( std::vector<HistogramBead>& bead, const std::vector<double>& gpoint, const std::vector<double>& args,
456 : const double& height, std::vector<double>& der ) const {
457 533600 : double val=height; std::vector<double> contr( args.size() ); Value* bw_arg=getPntrToArgument(bwargno);
458 533600 : if( bw_arg->getRank()<2 ) {
459 1067200 : for(unsigned j=0; j<args.size(); ++j) {
460 533600 : bead[j].set( gpoint[j], gpoint[j]+gridobject.getGridSpacing()[j], 1/sqrt(bw_arg->get(j)) );
461 533600 : contr[j] = bead[j].calculateWithCutoff( args[j], der[j] );
462 533600 : val = val*contr[j];
463 : }
464 0 : } else plumed_error();
465 1067200 : for(unsigned j=0; j<args.size(); ++j) {
466 533600 : if( fabs(contr[j])>epsilon ) der[j] *= val / contr[j];
467 : }
468 533600 : return val;
469 : }
470 :
471 251539 : void KDE::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
472 : const unsigned& bufstart, std::vector<double>& buffer ) const {
473 : plumed_dbg_assert( valindex==0 );
474 251539 : if( numberOfKernels==1 ) {
475 21277 : unsigned istart = bufstart + (1+gridobject.getDimension())*code;
476 21277 : unsigned valout = getConstPntrToComponent(0)->getPositionInStream(); buffer[istart] += myvals.get( valout );
477 78047 : for(unsigned i=0; i<gridobject.getDimension(); ++i) buffer[istart+1+i] += myvals.getDerivative( valout, i );
478 46591 : return;
479 : }
480 230262 : std::vector<double> args( gridobject.getDimension() ); double height; retrieveArgumentsAndHeight( myvals, args, height );
481 230262 : if( !ignore_out_of_bounds && !gridobject.inbounds( args ) ) {
482 : // if( fabs(height)>epsilon ) warning("bounds are possibly set too small as hills with substantial heights are being ignored");
483 : return ;
484 : }
485 : // Add the kernel to the grid
486 : unsigned num_neigh; std::vector<unsigned> neighbors;
487 204948 : if( kerneltype!="DISCRETE" ) gridobject.getNeighbors( args, nneigh, num_neigh, neighbors );
488 204948 : std::vector<double> der( args.size() ), gpoint( args.size() );
489 204948 : if( fabs(height)>epsilon ) {
490 204948 : if( getName()=="KDE" ) {
491 197340 : if( kerneltype=="DISCRETE" ) {
492 31494 : std::vector<double> newargs( args.size() );
493 62988 : for(unsigned i=0; i<args.size(); ++i) newargs[i] = args[i] + 0.5*gridobject.getGridSpacing()[i];
494 31494 : plumed_assert( bufstart + gridobject.getIndex( newargs )*(1+args.size())<buffer.size() );
495 31494 : buffer[ bufstart + gridobject.getIndex( newargs )*(1+args.size()) ] += height;
496 165846 : } else if( kerneltype.find("bin")!=std::string::npos ) {
497 104400 : std::vector<HistogramBead> bead( args.size() ); setupHistogramBeads( bead );
498 522000 : for(unsigned i=0; i<num_neigh; ++i) {
499 417600 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
500 417600 : double val = evaluateBeadValue( bead, gpoint, args, height, der );
501 417600 : buffer[ bufstart + neighbors[i]*(1+der.size()) ] += val;
502 835200 : for(unsigned j=0; j<der.size(); ++j) buffer[ bufstart + neighbors[i]*(1+der.size()) + 1 + j ] += val*der[j];
503 : }
504 : } else {
505 37979124 : for(unsigned i=0; i<num_neigh; ++i) {
506 37917678 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
507 37917678 : buffer[ bufstart + neighbors[i]*(1+der.size()) ] += evaluateKernel( gpoint, args, height, der );
508 150985777 : for(unsigned j=0; j<der.size(); ++j) buffer[ bufstart + neighbors[i]*(1+der.size()) + 1 + j ] += der[j];
509 : }
510 : }
511 : } else {
512 828405 : for(unsigned i=0; i<num_neigh; ++i) {
513 820797 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
514 3283188 : double dot=0; for(unsigned j=0; j<gpoint.size(); ++j) dot += args[j]*gpoint[j];
515 820797 : double newval = height*von_misses_norm*exp( von_misses_concentration*dot );
516 820797 : buffer[ bufstart + neighbors[i]*(1+gpoint.size()) ] += newval;
517 3283188 : for(unsigned j=0; j<gpoint.size(); ++j) buffer[ bufstart + neighbors[i]*(1+gpoint.size()) + 1 + j ] += von_misses_concentration*newval*gpoint[j];
518 : }
519 : }
520 : }
521 : }
522 :
523 610 : void KDE::updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const {
524 610 : if( !myval->forcesWereAdded() ) return ;
525 610 : if( numberOfKernels==1 ) plumed_error();
526 :
527 610 : int flag=1;
528 290135 : for(unsigned i=0; i<numberOfKernels; ++i) {
529 289525 : if( checkTaskStatus( i, flag ) ) force_tasks.push_back(i);
530 : }
531 : }
532 :
533 58028 : void KDE::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
534 58028 : if( numberOfKernels==1 ) {
535 0 : plumed_error();
536 : return;
537 : }
538 58028 : double height; std::vector<double> args( gridobject.getDimension() );
539 58028 : retrieveArgumentsAndHeight( myvals, args, height );
540 : unsigned num_neigh; std::vector<unsigned> neighbors;
541 58028 : gridobject.getNeighbors( args, nneigh, num_neigh, neighbors );
542 58028 : std::vector<double> der( args.size() ), gpoint( args.size() );
543 129700 : unsigned hforce_start = 0; for(unsigned j=0; j<der.size(); ++j) hforce_start += getPntrToArgument(j)->getNumberOfStoredValues();
544 58028 : if( fabs(height)>epsilon ) {
545 58028 : if( getName()=="KDE" ) {
546 51226 : if( kerneltype.find("bin")!=std::string::npos ) {
547 29000 : std::vector<HistogramBead> bead( args.size() ); setupHistogramBeads( bead );
548 145000 : for(unsigned i=0; i<num_neigh; ++i) {
549 116000 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
550 116000 : double val = evaluateBeadValue( bead, gpoint, args, height, der ); double fforce = getConstPntrToComponent(0)->getForce( neighbors[i] );
551 116000 : if( hasheight && getPntrToArgument(args.size())->getRank()==0 ) forces[ hforce_start ] += val*fforce / height;
552 116000 : else if( hasheight ) forces[ hforce_start + getPntrToArgument(args.size())->getIndexInStore(itask) ] += val*fforce / height;
553 232000 : unsigned n=0; for(unsigned j=0; j<der.size(); ++j) { forces[n + getPntrToArgument(j)->getIndexInStore(itask)] += der[j]*fforce; n += getPntrToArgument(j)->getNumberOfStoredValues(); }
554 : }
555 : } else {
556 202954 : for(unsigned i=0; i<num_neigh; ++i) {
557 180728 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
558 180728 : double val = evaluateKernel( gpoint, args, height, der ), fforce = getConstPntrToComponent(0)->getForce( neighbors[i] );
559 180728 : if( hasheight && getPntrToArgument(args.size())->getRank()==0 ) forces[ hforce_start ] += val*fforce / height;
560 179910 : else if( hasheight ) forces[ hforce_start + getPntrToArgument(args.size())->getIndexInStore(itask) ] += val*fforce / height;
561 364382 : unsigned n=0; for(unsigned j=0; j<der.size(); ++j) { forces[n + getPntrToArgument(j)->getIndexInStore(itask)] += -der[j]*fforce; n += getPntrToArgument(j)->getNumberOfStoredValues(); }
562 : }
563 : }
564 : } else {
565 687002 : for(unsigned i=0; i<num_neigh; ++i) {
566 680200 : gridobject.getGridPointCoordinates( neighbors[i], gpoint );
567 2720800 : double dot=0; for(unsigned j=0; j<gpoint.size(); ++j) dot += args[j]*gpoint[j];
568 680200 : double fforce = myval->getForce( neighbors[i] ); double newval = height*von_misses_norm*exp( von_misses_concentration*dot );
569 680200 : if( hasheight && getPntrToArgument(args.size())->getRank()==0 ) forces[ hforce_start ] += newval*fforce / height;
570 680200 : else if( hasheight ) forces[ hforce_start + getPntrToArgument(args.size())->getIndexInStore(itask) ] += newval*fforce / height;
571 2720800 : unsigned n=0; for(unsigned j=0; j<gpoint.size(); ++j) { forces[n + getPntrToArgument(j)->getIndexInStore(itask)] += von_misses_concentration*newval*gpoint[j]*fforce; n += getPntrToArgument(j)->getNumberOfStoredValues(); }
572 : }
573 : }
574 : }
575 : }
576 :
577 : }
578 : }
|
