Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "core/ActionShortcut.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/ActionWithValue.h"
27 : #include "tools/IFile.h"
28 :
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace refdist {
33 :
34 : //+PLUMEDOC FUNCTION KERNEL
35 : /*
36 : Transform a set of input coordinates using a kernel function
37 :
38 : This action takes a vector of arguments in input, $s$, the square of the [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) or
39 : [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) between the instataneous values of these arguments
40 : and a set of reference values for them is then computed. If this squared distance is $x$ then the final quantity that is returned
41 : by this function is:
42 :
43 : $$
44 : f = w K(x)
45 : $$
46 :
47 : where $w$ is a scalar that the user specifies using the `WEIGHT` keyword and $K$ can be a function of $x$ that the user
48 : specifies using the `FUNC` keyword or one of the Kernel functions options from the following table that are available internally:
49 :
50 : | Instruction | Function |
51 : |:-----------:|:---------|
52 : | FUNC=gaussian | $\exp(-x/2)$ |
53 : | FUNC=von-misses | $\exp(-x/2)$ |
54 : | FUNC=triangular | $1-\sqrt{x} \quad \textrm{if} \quad x<1 \quad \textrm{otherwise} \quad 0$ |
55 :
56 : The `von-misses` options here is used if the input arguments have a periodic domain. This instruction changes the way the
57 : [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) is computed so that the method that is appropriate for using with periodic
58 : variables is employed in place of the default method.
59 :
60 : ## Examples
61 :
62 : The following example demonstrates how this action can be used to evaluate a Kernel that is a function of one argument:
63 :
64 : ```plumed
65 : d: DISTANCE ATOMS=1,2
66 : k: KERNEL ARG=d TYPE=gaussian CENTER=1 SIGMA=0.1
67 : ```
68 :
69 : The [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) between the instantaneous distance and the point 1 is evaluated here.
70 : The metric vector that is used to evaluate this distance is taking the reciprocal of the square of the input SIGMA value. Lastly, the height of the Gaussian, $w$,
71 : in the expression above is set equal to one. If you would like the total volume of the Gaussian to be equal to one you use the `NORMALIZED` keyword as shown here:
72 :
73 : ```plumed
74 : d: DISTANCE ATOMS=1,2
75 : k: KERNEL ARG=d TYPE=gaussian CENTER=1 SIGMA=0.1 NORMALIZED
76 : ```
77 :
78 : If your Kernel is a function of two arguments you can use the [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) to evaluate the
79 : distance as is done in the following example:
80 :
81 : ```plumed
82 : d1: DISTANCE ATOMS=1,2
83 : d2: DISTANCE ATOMS=3,4
84 : k: KERNEL ARG=d1,d2 TYPE=gaussian CENTER=1,1 SIGMA=0.1,0.1
85 : ```
86 :
87 : or you can use the [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) as is done here:
88 :
89 : ```plumed
90 : phi: TORSION ATOMS=1,2,3,4
91 : psi: TORSION ATOMS=5,6,7,8
92 : k: KERNEL ...
93 : ARG=phi,psi TYPE=von-misses
94 : CENTER=-1.09648066E+0000,-7.17867907E-0001
95 : COVAR=1.40523052E-0001,-1.05385552E-0001,-1.05385552E-0001,1.63290557E-0001
96 : ...
97 : ```
98 :
99 : To switch to using the [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) you use `COVAR` instead of `SIGMA` and specify the covariance matrix of the kernel.
100 : The metric that is used when evaluating the distance is the inverse of the input covariance matrix.
101 :
102 : Notice that you specify the Kernel function in a separate PLUMED input file by using an input like the one shown below:
103 :
104 : ```plumed
105 : #SETTINGS INPUTFILES=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp
106 : phi: TORSION ATOMS=1,2,3,4
107 : psi: TORSION ATOMS=5,6,7,8
108 : k: KERNEL ARG=phi,psi REFERENCE=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp NUMBER=3
109 : ```
110 :
111 : This command computes the same Kernel that was computed in the previous input. The keyword `REFERENCE` specifies that the parameters are to be read
112 : from the file and the keyword `NUMBER` indicates that the parameters are found on the third line of that file.
113 :
114 : Lastly, note that you can use vectors in the input to this shortcut as shown here:
115 :
116 : ```plumed
117 : #SETTINGS INPUTFILES=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp
118 : phi: TORSION ATOMS1=1,2,3,4 ATOMS2=9,10,11,12 ATOMS3=17,18,19,20
119 : psi: TORSION ATOMS1=5,6,7,8 ATOMS2=13,14,15,16 ATOMS3=21,22,23,24
120 : k: KERNEL ARG=phi,psi REFERENCE=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp NUMBER=3
121 : ```
122 :
123 : The output `k` here is a vector with three components. The first component is the kernel evaluated with the torsion involving atoms 1, 2, 3 and 4 and the
124 : torsion involving atoms 5, 6, 7 and 8. The second component is the kernel evaluated with the torsion involving atoms 9, 10, 11 and 12 and the
125 : torsion involving atoms 13, 14, 15 and 16. The third component is the kernel evaluated with the torsion involving atoms 17, 18, 19 and 20 and the
126 : torsion involving atoms 21, 22, 23 and 24.
127 :
128 : */
129 : //+ENDPLUMEDOC
130 :
131 :
132 : class Kernel : public ActionShortcut {
133 : public:
134 : static std::string fixArgumentDot( const std::string& argin );
135 : explicit Kernel(const ActionOptions&);
136 : static void registerKeywords(Keywords& keys);
137 : };
138 :
139 :
140 : PLUMED_REGISTER_ACTION(Kernel,"KERNEL")
141 :
142 20 : void Kernel::registerKeywords(Keywords& keys) {
143 20 : ActionShortcut::registerKeywords( keys );
144 40 : keys.addInputKeyword("numbered","ARG","scalar/vector","the arguments that should be used as input to this method");
145 20 : keys.add("compulsory","TYPE","gaussian","the type of kernel to use");
146 20 : keys.add("compulsory","CENTER","the position of the center of the kernel");
147 20 : keys.add("optional","SIGMA","square root of variance of the cluster");
148 20 : keys.add("compulsory","COVAR","the covariance of the kernel");
149 20 : keys.add("compulsory","WEIGHT","1.0","the weight to multiply this kernel function by");
150 20 : keys.add("optional","REFERENCE","the file from which to read the kernel parameters");
151 20 : keys.add("compulsory","NUMBER","1","if there are multiple sets of kernel parameters in the input file which set of kernel parameters would you like to read in here");
152 20 : keys.addFlag("NORMALIZED",false,"would you like the kernel function to be normalized");
153 40 : keys.setValueDescription("scalar/vector","the value of the kernel evaluated at the argument values");
154 20 : keys.needsAction("CONSTANT");
155 20 : keys.needsAction("CUSTOM");
156 20 : keys.needsAction("NORMALIZED_EUCLIDEAN_DISTANCE");
157 20 : keys.needsAction("PRODUCT");
158 20 : keys.needsAction("INVERT_MATRIX");
159 20 : keys.needsAction("MAHALANOBIS_DISTANCE");
160 20 : keys.needsAction("DIAGONALIZE");
161 20 : keys.needsAction("CONCATENATE");
162 20 : keys.needsAction("DETERMINANT");
163 20 : keys.needsAction("BESSEL");
164 20 : }
165 :
166 32 : std::string Kernel::fixArgumentDot( const std::string& argin ) {
167 32 : std::string argout = argin;
168 32 : std::size_t dot=argin.find(".");
169 32 : if( dot!=std::string::npos ) {
170 0 : argout = argin.substr(0,dot) + "_" + argin.substr(dot+1);
171 : }
172 32 : return argout;
173 : }
174 :
175 9 : Kernel::Kernel(const ActionOptions&ao):
176 : Action(ao),
177 9 : ActionShortcut(ao) {
178 : // Read in the arguments
179 : std::vector<std::string> argnames;
180 18 : parseVector("ARG",argnames);
181 9 : if( argnames.size()==0 ) {
182 0 : error("no arguments were specified");
183 : }
184 : // Now sort out the parameters
185 : double weight;
186 : std::string fname;
187 18 : parse("REFERENCE",fname);
188 : bool usemahalanobis=false;
189 9 : if( fname.length()>0 ) {
190 9 : IFile ifile;
191 9 : ifile.open(fname);
192 9 : ifile.allowIgnoredFields();
193 : unsigned number;
194 9 : parse("NUMBER",number);
195 : bool readline=false;
196 : // Create actions to hold the position of the center
197 31 : for(unsigned line=0; line<number; ++line) {
198 90 : for(unsigned i=0; i<argnames.size(); ++i) {
199 : std::string val;
200 59 : ifile.scanField(argnames[i], val);
201 59 : if( line==number-1 ) {
202 32 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref: CONSTANT VALUES=" + val );
203 : }
204 : }
205 62 : if( ifile.FieldExist("sigma_" + argnames[0]) ) {
206 : std::string varstr;
207 0 : for(unsigned i=0; i<argnames.size(); ++i) {
208 : std::string val;
209 0 : ifile.scanField("sigma_" + argnames[i], val);
210 0 : if( i==0 ) {
211 : varstr = val;
212 : } else {
213 0 : varstr += "," + val;
214 : }
215 : }
216 0 : if( line==number-1 ) {
217 0 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + varstr );
218 : }
219 : } else {
220 : std::string varstr, nvals;
221 31 : Tools::convert( argnames.size(), nvals );
222 31 : usemahalanobis=(argnames.size()>1);
223 90 : for(unsigned i=0; i<argnames.size(); ++i) {
224 174 : for(unsigned j=0; j<argnames.size(); j++) {
225 : std::string val;
226 230 : ifile.scanField("sigma_" +argnames[i] + "_" + argnames[j], val );
227 115 : if(i==0 && j==0 ) {
228 : varstr = val;
229 : } else {
230 168 : varstr += "," + val;
231 : }
232 : }
233 : }
234 31 : if( line==number-1 ) {
235 9 : if( !usemahalanobis ) {
236 4 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + varstr );
237 : } else {
238 14 : readInputLine( getShortcutLabel() + "_cov: CONSTANT NCOLS=" + nvals + " NROWS=" + nvals + " VALUES=" + varstr );
239 : }
240 : }
241 : }
242 31 : if( line==number-1 ) {
243 : readline=true;
244 : break;
245 : }
246 22 : ifile.scanField();
247 : }
248 9 : if( !readline ) {
249 0 : error("could not read reference configuration");
250 : }
251 9 : ifile.scanField();
252 9 : ifile.close();
253 9 : } else {
254 : // Create actions to hold the position of the center
255 0 : std::vector<std::string> center(argnames.size());
256 0 : parseVector("CENTER",center);
257 0 : for(unsigned i=0; i<argnames.size(); ++i) {
258 0 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref: CONSTANT VALUES=" + center[i] );
259 : }
260 : std::vector<std::string> sig;
261 0 : parseVector("SIGMA",sig);
262 0 : if( sig.size()==0 ) {
263 : // Create actions to hold the covariance
264 : std::string cov;
265 0 : parse("COVAR",cov);
266 0 : usemahalanobis=(argnames.size()>1);
267 0 : if( !usemahalanobis ) {
268 0 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + cov );
269 : } else {
270 : std::string nvals;
271 0 : Tools::convert( argnames.size(), nvals );
272 0 : readInputLine( getShortcutLabel() + "_cov: CONSTANT NCOLS=" + nvals + " NROWS=" + nvals + " VALUES=" + cov );
273 : }
274 0 : } else if( sig.size()==argnames.size() ) {
275 : // And actions to hold the standard deviation
276 0 : std::string valstr = sig[0];
277 0 : for(unsigned i=1; i<sig.size(); ++i) {
278 0 : valstr += "," + sig[i];
279 : }
280 0 : readInputLine( getShortcutLabel() + "_sigma: CONSTANT VALUES=" + valstr );
281 0 : readInputLine( getShortcutLabel() + "_var: CUSTOM ARG=" + getShortcutLabel() + "_sigma FUNC=x*x PERIODIC=NO");
282 : } else {
283 0 : error("sigma has wrong length");
284 : }
285 0 : }
286 :
287 : // Create the reference point and arguments
288 : std::string refpoint, argstr;
289 25 : for(unsigned i=0; i<argnames.size(); ++i) {
290 16 : if( i==0 ) {
291 : argstr = argnames[0];
292 18 : refpoint = getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref";
293 : } else {
294 14 : argstr += "," + argnames[1];
295 14 : refpoint += "," + getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref";
296 : }
297 : }
298 :
299 : // Get the information on the kernel type
300 : std::string func_str, ktype;
301 18 : parse("TYPE",ktype);
302 16 : if( ktype=="gaussian" || ktype=="von-misses" ) {
303 : func_str = "exp(-x/2)";
304 0 : } else if( ktype=="triangular" ) {
305 : func_str = "step(1.-sqrt(x))*(1.-sqrt(x))";
306 : } else {
307 : func_str = ktype;
308 : }
309 9 : std::string vm_str="";
310 9 : if( ktype=="von-misses" ) {
311 : vm_str=" VON_MISSES";
312 : }
313 :
314 9 : unsigned nvals = argnames.size();
315 : bool norm;
316 9 : parseFlag("NORMALIZED",norm);
317 9 : if( !usemahalanobis ) {
318 : // Invert the variance
319 4 : readInputLine( getShortcutLabel() + "_icov: CUSTOM ARG=" + getShortcutLabel() + "_var FUNC=1/x PERIODIC=NO");
320 : // Compute the distance between the center of the basin and the current configuration
321 4 : readInputLine( getShortcutLabel() + "_dist_2: NORMALIZED_EUCLIDEAN_DISTANCE SQUARED" + vm_str +" ARG1=" + argstr + " ARG2=" + refpoint + " METRIC=" + getShortcutLabel() + "_icov");
322 : // And compute a determinent for the input covariance matrix if it is required
323 2 : if( norm ) {
324 2 : if( ktype=="von-misses" ) {
325 0 : readInputLine( getShortcutLabel() + "_vec: CUSTOM ARG=" + getShortcutLabel() + "_icov FUNC=x PERIODIC=NO" );
326 : } else {
327 4 : readInputLine( getShortcutLabel() + "_det: PRODUCT ARG=" + getShortcutLabel() + "_var");
328 : }
329 : }
330 : } else {
331 : // Invert the input covariance matrix
332 14 : readInputLine( getShortcutLabel() + "_icov: INVERT_MATRIX ARG=" + getShortcutLabel() + "_cov" );
333 : // Compute the distance between the center of the basin and the current configuration
334 14 : readInputLine( getShortcutLabel() + "_dist_2: MAHALANOBIS_DISTANCE SQUARED ARG1=" + argstr + " ARG2=" + refpoint + " METRIC=" + getShortcutLabel() + "_icov " + vm_str );
335 : // And compute a determinent for the input covariance matrix if it is required
336 7 : if( norm ) {
337 7 : if( ktype=="von-misses" ) {
338 14 : readInputLine( getShortcutLabel() + "_det: DIAGONALIZE ARG=" + getShortcutLabel() + "_cov VECTORS=all" );
339 7 : std::string num, argnames= getShortcutLabel() + "_det.vals-1";
340 14 : for(unsigned i=1; i<nvals; ++i) {
341 7 : Tools::convert( i+1, num );
342 14 : argnames += "," + getShortcutLabel() + "_det.vals-" + num;
343 : }
344 14 : readInputLine( getShortcutLabel() + "_comp: CONCATENATE ARG=" + argnames );
345 14 : readInputLine( getShortcutLabel() + "_vec: CUSTOM ARG=" + getShortcutLabel() + "_comp FUNC=1/x PERIODIC=NO");
346 : } else {
347 0 : readInputLine( getShortcutLabel() + "_det: DETERMINANT ARG=" + getShortcutLabel() + "_cov");
348 : }
349 : }
350 : }
351 :
352 : // Compute the Gaussian
353 : std::string wstr;
354 9 : parse("WEIGHT",wstr);
355 9 : if( norm ) {
356 9 : if( ktype=="gaussian" ) {
357 : std::string pstr;
358 2 : Tools::convert( sqrt(pow(2*pi,nvals)), pstr );
359 4 : readInputLine( getShortcutLabel() + "_vol: CUSTOM ARG=" + getShortcutLabel() + "_det FUNC=(sqrt(x)*" + pstr + ") PERIODIC=NO");
360 7 : } else if( ktype=="von-misses" ) {
361 : std::string wstr, min, max;
362 7 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_dist_2_diff" );
363 7 : plumed_assert( av );
364 7 : if( !av->copyOutput(0)->isPeriodic() ) {
365 0 : error("VON_MISSES only works with periodic variables");
366 : }
367 7 : av->copyOutput(0)->getDomain(min,max);
368 14 : readInputLine( getShortcutLabel() + "_bes: BESSEL ORDER=0 ARG=" + getShortcutLabel() + "_vec");
369 14 : readInputLine( getShortcutLabel() + "_cc: CUSTOM ARG=" + getShortcutLabel() + "_bes FUNC=("+max+"-"+min+")*x PERIODIC=NO");
370 14 : readInputLine( getShortcutLabel() + "_vol: PRODUCT ARG=" + getShortcutLabel() + "_cc");
371 : } else {
372 0 : error("only gaussian and von-misses kernels are normalizable");
373 : }
374 : // And the (suitably normalized) kernel
375 18 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_dist_2," + getShortcutLabel() + "_vol FUNC=" + wstr + "*exp(-x/2)/y PERIODIC=NO");
376 : } else {
377 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_dist_2 FUNC=" + wstr + "*" + func_str + " PERIODIC=NO");
378 : }
379 9 : checkRead();
380 :
381 9 : }
382 :
383 : }
384 : }
385 :
386 :
|
