Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 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_KDEGridTools_h
23 : #define __PLUMED_gridtools_KDEGridTools_h
24 :
25 : #include "KDE.h"
26 :
27 : namespace PLMD {
28 : namespace gridtools {
29 :
30 : template <class K, class P>
31 : class KDEGridTools {
32 : public:
33 : std::vector<double> gspacing;
34 : std::vector<std::size_t> nbin;
35 : std::vector<std::string> gmin, gmax;
36 : static void registerKeywords( Keywords& keys );
37 : static void readHeightKeyword( bool canusevol, std::size_t nargs, const std::vector<std::string>& bw, ActionWithArguments* action );
38 : static void readBandwidth( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw );
39 : static void readBandwidthKeyword( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw, std::vector<Value*>& bwargs );
40 : static void readBandwidthAndHeight( const P& params, ActionWithArguments* action );
41 : static void convertHeightsToVolumes( const std::size_t& nargs, const std::vector<std::string>& bw, const std::string& volstr, ActionWithArguments* action );
42 : static void readGridParameters( KDEGridTools<K,P>& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape );
43 : static void setupGridBounds( KDEGridTools<K,P>& g, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval );
44 : static void getDiscreteSupport( const KDEGridTools<K,P>& g, P& p, const View<const double>& shape, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject );
45 : static void getNeighbors( const P& p, View<double> at, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors );
46 : };
47 :
48 : template <class K, class P>
49 397 : void KDEGridTools<K,P>::registerKeywords( Keywords& keys ) {
50 397 : keys.add("optional","BANDWIDTH","the bandwidths for kernel density esimtation");
51 397 : keys.add("optional","VOLUMES","this keyword take the label of an action that calculates a vector of values. The elements of this vector "
52 : "divided by the volume of the Gaussian are used as weights for the Gaussians");
53 397 : keys.add("optional","HEIGHTS","this keyword takes the label of an action that calculates a vector of values. The elements of this vector "
54 : "are used as weights for the Gaussians.");
55 397 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
56 397 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
57 397 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
58 397 : keys.add("optional","GRID_BIN","the number of bins for the grid");
59 397 : }
60 :
61 : template <class K, class P>
62 80 : void KDEGridTools<K,P>::readHeightKeyword( bool canusevol, std::size_t nargs, const std::vector<std::string>& bw, ActionWithArguments* action ) {
63 : std::string weight_str;
64 160 : action->parse("HEIGHTS",weight_str);
65 : std::string str_nvals;
66 80 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
67 : std::string nr, nc;
68 4 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
69 4 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
70 8 : str_nvals = nr + "," + nc;
71 : } else {
72 76 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), str_nvals );
73 : }
74 80 : if( weight_str.length()>0 ) {
75 38 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( weight_str, action, "_heights", true );
76 61 : } else if( canusevol ) {
77 40 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_volumes: ONES SIZE=" + str_nvals ), false );
78 80 : KDEGridTools<K,P>::convertHeightsToVolumes(nargs,bw,action->getLabel() + "_volumes",action);
79 : } else {
80 21 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: ONES SIZE=" + str_nvals ), false );
81 42 : KDEHelper<K,P,KDEGridTools<K,P>>::addArgument( action->getLabel() + "_heights", action );
82 : }
83 82 : }
84 :
85 : template<class K, class P>
86 84 : void KDEGridTools<K, P>::readBandwidthKeyword( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw, std::vector<Value*>& bwargs ) {
87 84 : action->parseVector("BANDWIDTH",bw);
88 84 : if( nargs>1 && bw.size()==1 ) {
89 0 : ActionWithArguments::interpretArgumentList( bw, action->plumed.getActionSet(), action, bwargs );
90 0 : if( bwargs.size()!=1 ) {
91 0 : action->error("invalid bandwidth found");
92 : }
93 : // Create a vector of ones with the right size
94 : std::string nvals;
95 0 : if( (action->getPntrToArgument(0))->getRank()==2 ) {
96 : std::string nr, nc;
97 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[0], nr );
98 0 : Tools::convert( (action->getPntrToArgument(0))->getShape()[1], nc );
99 0 : nvals = nr + "," + nc;
100 : } else {
101 0 : Tools::convert( (action->getPntrToArgument(0))->getNumberOfValues(), nvals );
102 : }
103 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_bwones: ONES SIZE=" + nvals ), false );
104 : }
105 84 : }
106 :
107 : template <class K, class P>
108 83 : void KDEGridTools<K,P>::readBandwidth( std::size_t nargs, ActionWithArguments* action, std::vector<std::string>& bw ) {
109 83 : plumed_assert( typeid(K)==typeid(DiagonalKernelParams) );
110 83 : std::vector<Value*> bwargs;
111 83 : readBandwidthKeyword( nargs, action, bw, bwargs );
112 83 : if( nargs>1 && bw.size()==1 ) {
113 0 : if( bwargs[0]->getRank()!=1 || bwargs[0]->getNumberOfValues()!=nargs ) {
114 0 : action->error("invalid input for bandwidth parameter");
115 : }
116 : std::string str_i;
117 0 : bw.resize( nargs );
118 0 : for(unsigned i=0; i<nargs; ++i) {
119 0 : Tools::convert( i+1, str_i );
120 0 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_scalar_bw" + str_i + ": SELECT_COMPONENTS ARG=" + bwargs[0]->getName() + " COMPONENTS=" + str_i ), false );
121 0 : bw[i] = action->getLabel() + "_bw" + str_i;
122 0 : action->plumed.readInputWords( Tools::getWords( bw[i] + ": CUSTOM ARG=" + action->getLabel() + "_bwones," + action->getLabel() + "_scalar_bw" + str_i + " FUNC=x*y PERIODIC=NO"), false );
123 : }
124 83 : } else if( bw.size()==nargs ) {
125 : double bwval;
126 : std::string str_i;
127 197 : for(unsigned i=0; i<nargs; ++i) {
128 115 : Tools::convert( i+1, str_i );
129 115 : if( Tools::convertNoexcept( bw[i], bwval ) && bwval<0 ) {
130 1 : action->error("BANDWIDTH should not be negative");
131 : }
132 114 : if( Tools::convertNoexcept( bw[i], bwval ) && fabs(bwval)<epsilon ) {
133 4 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( bw[i], action, "_bwz" + str_i, true );
134 : } else {
135 225 : KDEHelper<K,P,KDEGridTools<K,P>>::readKernelParameters( bw[i], action, "_bw" + str_i, true );
136 : }
137 : }
138 : } else {
139 1 : action->error("wrong number of arguments specified in input to bandwidth parameter");
140 : }
141 82 : }
142 :
143 : template <class K, class P>
144 58 : void KDEGridTools<K, P>::convertHeightsToVolumes( const std::size_t& nargs, const std::vector<std::string>& bw, const std::string& volstr, ActionWithArguments* action ) {
145 58 : if( bw.size()==nargs ) {
146 58 : unsigned nonzeroargs=0;
147 140 : for(unsigned i=0; i<nargs; ++i) {
148 82 : if( bw[i].find("_bwz")==std::string::npos ) {
149 80 : nonzeroargs++;
150 : }
151 : }
152 : std::string str_i, nargs_str;
153 58 : Tools::convert( nonzeroargs, nargs_str );
154 116 : std::string varstr = "VAR=h", funcstr = "FUNC=h/(sqrt((2*pi)^" + nargs_str + ")", argstr = "ARG=" + volstr;
155 140 : for(unsigned i=0; i<nargs; ++i) {
156 82 : if( bw[i].find("_bwz")!=std::string::npos ) {
157 2 : continue;
158 : }
159 80 : Tools::convert( i+1, str_i );
160 80 : varstr += ",b" + str_i;
161 160 : funcstr += "*b" + str_i;
162 160 : argstr += "," + bw[i];
163 : }
164 : funcstr += ")";
165 58 : if( (action->getPntrToArgument(0))->getNumberOfValues()==1 ) {
166 152 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO " + argstr + " " + varstr + " " + funcstr), false );
167 : } else {
168 100 : action->plumed.readInputWords( Tools::getWords(action->getLabel() + "_heights: CUSTOM PERIODIC=NO " + argstr + " " + varstr + " " + funcstr + " MASK=" + volstr), false );
169 : }
170 116 : KDEHelper<K,P,KDEGridTools<K,P>>::addArgument( action->getLabel() + "_heights", action );
171 0 : } else if( bw.size()==nargs*nargs ) {
172 0 : action->error("have not implemented normalization parameter for non-diagonal kernels");
173 : }
174 :
175 58 : }
176 :
177 : template <class K, class P>
178 98 : void KDEGridTools<K,P>::readGridParameters( KDEGridTools<K,P>& g, ActionWithArguments* action, GridCoordinatesObject& gridobject, std::vector<std::size_t>& shape ) {
179 98 : g.gmin.resize( shape.size() );
180 98 : g.gmax.resize( shape.size() );
181 98 : action->parseVector("GRID_MIN",g.gmin);
182 98 : action->parseVector("GRID_MAX",g.gmax);
183 230 : for(unsigned i=0; i<g.gmin.size(); ++i) {
184 132 : if( g.gmin[i]=="auto" ) {
185 52 : action->log.printf(" for %dth coordinate min and max are set automatically \n", (i+1) );
186 52 : if( g.gmax[i]!="auto" ) {
187 0 : action->error("if gmin is set automatically gmax must also be set automatically");
188 : }
189 : Value* myarg = action->getPntrToArgument(i);
190 52 : if( myarg->isPeriodic() ) {
191 0 : if( g.gmin[i]=="auto" ) {
192 0 : myarg->getDomain( g.gmin[i], g.gmax[i] );
193 : } else {
194 : std::string str_min, str_max;
195 0 : myarg->getDomain( str_min, str_max );
196 0 : if( str_min!=g.gmin[i] || str_max!=g.gmax[i] ) {
197 0 : action->error("all periodic arguments should have the same domain");
198 : }
199 : }
200 52 : } else if( myarg->getName().find(".")!=std::string::npos ) {
201 52 : std::size_t dot = myarg->getName().find_first_of(".");
202 52 : std::string name = myarg->getName().substr(dot+1);
203 52 : if( name!="x" && name!="y" && name!="z" ) {
204 0 : action->error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
205 : }
206 : } else {
207 0 : action->error("cannot set GRID_MIN and GRID_MAX automatically if input argument is not component of distance");
208 : }
209 : } else {
210 80 : action->log.printf(" for %dth coordinate min is set to %s and max is set to %s \n", (i+1), g.gmin[i].c_str(), g.gmax[i].c_str() );
211 : }
212 : }
213 :
214 98 : action->parseVector("GRID_BIN",g.nbin);
215 196 : action->parseVector("GRID_SPACING",g.gspacing);
216 :
217 98 : if( g.nbin.size()!=shape.size() && g.gspacing.size()!=shape.size() ) {
218 0 : action->error("GRID_BIN or GRID_SPACING must be set");
219 : }
220 : // Create a value
221 98 : std::vector<bool> ipbc( shape.size() );
222 230 : for(unsigned i=0; i<shape.size(); ++i) {
223 132 : if( (action->getPntrToArgument( i ))->isPeriodic() || g.gmin[i]=="auto" ) {
224 : ipbc[i]=true;
225 67 : if( g.nbin.size()==shape.size() ) {
226 67 : shape[i] = g.nbin[i];
227 : }
228 : } else {
229 : ipbc[i]=false;
230 65 : if( g.nbin.size()==shape.size() ) {
231 64 : shape[i] = g.nbin[i]+1;
232 : }
233 : }
234 : }
235 196 : gridobject.setup( "flat", ipbc, 0, 0.0 );
236 98 : }
237 :
238 : template <class K, class P>
239 94 : void KDEGridTools<K,P>::setupGridBounds( KDEGridTools<K,P>& g, const Tensor& box, GridCoordinatesObject& gridobject, const std::vector<Value*>& args, Value* myval ) {
240 218 : for(unsigned i=0; i<gridobject.getDimension(); ++i) {
241 124 : if( g.gmin[i]=="auto" ) {
242 : double lcoord, ucoord;
243 46 : std::size_t dot = args[i]->getName().find_first_of(".");
244 46 : std::string name = args[i]->getName().substr(dot+1);
245 46 : if( name=="x" ) {
246 24 : lcoord=-0.5*box(0,0);
247 24 : ucoord=0.5*box(0,0);
248 22 : } else if( name=="y" ) {
249 12 : lcoord=-0.5*box(1,1);
250 12 : ucoord=0.5*box(1,1);
251 10 : } else if( name=="z" ) {
252 10 : lcoord=-0.5*box(2,2);
253 10 : ucoord=0.5*box(2,2);
254 : } else {
255 0 : plumed_error();
256 : }
257 : // And convert to strings for bin and bmax
258 46 : Tools::convert( lcoord, g.gmin[i] );
259 46 : Tools::convert( ucoord, g.gmax[i] );
260 : }
261 : }
262 : // And setup the grid object
263 94 : gridobject.setBounds( g.gmin, g.gmax, g.nbin, g.gspacing );
264 91 : myval->setShape( gridobject.getNbin(true) );
265 91 : }
266 :
267 : template <> inline
268 : void KDEGridTools<DiagonalKernelParams,DiscreteKernel>::getDiscreteSupport( const KDEGridTools<DiagonalKernelParams,DiscreteKernel>& g, DiscreteKernel& p, const View<const double>& shape, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
269 : return;
270 : }
271 :
272 : template <class K, class P>
273 77 : void KDEGridTools<K, P>::getDiscreteSupport( const KDEGridTools<K,P>& g, P& p, const View<const double>& shape, std::vector<unsigned>& nneigh, GridCoordinatesObject& gridobject ) {
274 77 : std::size_t ng = gridobject.getDimension();
275 77 : plumed_assert( nneigh.size()==ng );
276 77 : std::vector<double> support( ng );
277 77 : P::getSupport( p, shape, support );
278 183 : for(unsigned i=0; i<ng; ++i) {
279 106 : nneigh[i] = static_cast<unsigned>( ceil( support[i]/gridobject.getGridSpacing()[i] ));
280 : }
281 77 : }
282 :
283 : template <> inline
284 32295 : void KDEGridTools<DiagonalKernelParams,DiscreteKernel>::getNeighbors( const DiscreteKernel& p, View<double> at, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
285 32295 : num_neighbors=1;
286 32295 : neighbors.resize(1);
287 65374 : for(unsigned i=0; i<at.size(); ++i) {
288 33079 : at[i] += 0.5*gridobject.getGridSpacing()[i];
289 : }
290 32295 : neighbors[0]=gridobject.getIndex( View<const double>(at.data(),at.size()) );
291 32294 : }
292 :
293 : template <class K, class P>
294 211512 : void KDEGridTools<K,P>::getNeighbors( const P& p, View<double> at, const GridCoordinatesObject& gridobject, const std::vector<unsigned>& nneigh, unsigned& num_neighbors, std::vector<unsigned>& neighbors ) {
295 211512 : gridobject.getNeighbors( View<const double>(at.data(),at.size()), nneigh, num_neighbors, neighbors );
296 211512 : }
297 :
298 : }
299 : }
300 : #endif
|
