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