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 22 : void Kernel::registerKeywords(Keywords& keys) {
143 22 : ActionShortcut::registerKeywords( keys );
144 44 : keys.addInputKeyword("numbered","ARG","scalar/vector","the arguments that should be used as input to this method");
145 22 : keys.add("compulsory","TYPE","gaussian","the type of kernel to use");
146 22 : keys.add("compulsory","CENTER","the position of the center of the kernel");
147 22 : keys.add("optional","SIGMA","square root of variance of the cluster");
148 22 : keys.add("compulsory","COVAR","the covariance of the kernel");
149 22 : keys.add("compulsory","WEIGHT","1.0","the weight to multiply this kernel function by");
150 22 : keys.add("optional","REFERENCE","the file from which to read the kernel parameters");
151 22 : 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 22 : keys.addFlag("NORMALIZED",false,"would you like the kernel function to be normalized");
153 44 : keys.setValueDescription("scalar/vector","the value of the kernel evaluated at the argument values");
154 22 : keys.needsAction("CONSTANT");
155 22 : keys.needsAction("CUSTOM");
156 22 : keys.needsAction("NORMALIZED_EUCLIDEAN_DISTANCE");
157 22 : keys.needsAction("PRODUCT");
158 22 : keys.needsAction("INVERT_MATRIX");
159 22 : keys.needsAction("MAHALANOBIS_DISTANCE");
160 22 : keys.needsAction("DIAGONALIZE");
161 22 : keys.needsAction("CONCATENATE");
162 22 : keys.needsAction("DETERMINANT");
163 22 : keys.needsAction("BESSEL");
164 22 : }
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 9 : 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 : std::string fname;
186 18 : parse("REFERENCE",fname);
187 : bool usemahalanobis=false;
188 9 : if( fname.length()>0 ) {
189 9 : IFile ifile;
190 9 : ifile.open(fname);
191 9 : ifile.allowIgnoredFields();
192 : unsigned number;
193 9 : parse("NUMBER",number);
194 : bool readline=false;
195 : // Create actions to hold the position of the center
196 31 : for(unsigned line=0; line<number; ++line) {
197 90 : for(unsigned i=0; i<argnames.size(); ++i) {
198 : std::string val;
199 59 : ifile.scanField(argnames[i], val);
200 59 : if( line==number-1 ) {
201 48 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref: CONSTANT VALUES=" + val );
202 : }
203 : }
204 62 : if( ifile.FieldExist("sigma_" + argnames[0]) ) {
205 : std::string varstr;
206 0 : for(unsigned i=0; i<argnames.size(); ++i) {
207 : std::string val;
208 0 : ifile.scanField("sigma_" + argnames[i], val);
209 0 : if( i==0 ) {
210 : varstr = val;
211 : } else {
212 0 : varstr += "," + val;
213 : }
214 : }
215 0 : if( line==number-1 ) {
216 0 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + varstr );
217 : }
218 : } else {
219 : std::string varstr, nvals;
220 31 : Tools::convert( argnames.size(), nvals );
221 31 : usemahalanobis=(argnames.size()>1);
222 90 : for(unsigned i=0; i<argnames.size(); ++i) {
223 174 : for(unsigned j=0; j<argnames.size(); j++) {
224 : std::string val;
225 230 : ifile.scanField("sigma_" +argnames[i] + "_" + argnames[j], val );
226 115 : if(i==0 && j==0 ) {
227 : varstr = val;
228 : } else {
229 84 : varstr += "," + val;
230 : }
231 : }
232 : }
233 31 : if( line==number-1 ) {
234 9 : if( !usemahalanobis ) {
235 4 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + varstr );
236 : } else {
237 28 : readInputLine( getShortcutLabel() + "_cov: CONSTANT NCOLS=" + nvals + " NROWS=" + nvals + " VALUES=" + varstr );
238 : }
239 : }
240 : }
241 31 : if( line==number-1 ) {
242 : readline=true;
243 : break;
244 : }
245 22 : ifile.scanField();
246 : }
247 9 : if( !readline ) {
248 0 : error("could not read reference configuration");
249 : }
250 9 : ifile.scanField();
251 9 : ifile.close();
252 9 : } else {
253 : // Create actions to hold the position of the center
254 0 : std::vector<std::string> center(argnames.size());
255 0 : parseVector("CENTER",center);
256 0 : for(unsigned i=0; i<argnames.size(); ++i) {
257 0 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref: CONSTANT VALUES=" + center[i] );
258 : }
259 0 : std::vector<std::string> sig;
260 0 : parseVector("SIGMA",sig);
261 0 : if( sig.size()==0 ) {
262 : // Create actions to hold the covariance
263 : std::string cov;
264 0 : parse("COVAR",cov);
265 0 : usemahalanobis=(argnames.size()>1);
266 0 : if( !usemahalanobis ) {
267 0 : readInputLine( getShortcutLabel() + "_var: CONSTANT VALUES=" + cov );
268 : } else {
269 : std::string nvals;
270 0 : Tools::convert( argnames.size(), nvals );
271 0 : readInputLine( getShortcutLabel() + "_cov: CONSTANT NCOLS=" + nvals + " NROWS=" + nvals + " VALUES=" + cov );
272 : }
273 0 : } else if( sig.size()==argnames.size() ) {
274 : // And actions to hold the standard deviation
275 0 : std::string valstr = sig[0];
276 0 : for(unsigned i=1; i<sig.size(); ++i) {
277 0 : valstr += "," + sig[i];
278 : }
279 0 : readInputLine( getShortcutLabel() + "_sigma: CONSTANT VALUES=" + valstr );
280 0 : readInputLine( getShortcutLabel() + "_var: CUSTOM ARG=" + getShortcutLabel() + "_sigma FUNC=x*x PERIODIC=NO");
281 : } else {
282 0 : error("sigma has wrong length");
283 : }
284 0 : }
285 :
286 : // Create the reference point and arguments
287 : std::string refpoint;
288 : std::string argstr;
289 25 : for(unsigned i=0; i<argnames.size(); ++i) {
290 16 : if( i==0 ) {
291 : argstr = argnames[0];
292 27 : refpoint = getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref";
293 : } else {
294 14 : argstr += "," + argnames[1];
295 21 : refpoint += "," + getShortcutLabel() + "_" + fixArgumentDot(argnames[i]) + "_ref";
296 : }
297 : }
298 :
299 : // Get the information on the kernel type
300 : std::string func_str;
301 : std::string ktype;
302 9 : parse("TYPE",ktype);
303 9 : if( ktype=="gaussian" || ktype=="von-misses" ) {
304 : func_str = "exp(-x/2)";
305 0 : } else if( ktype=="triangular" ) {
306 : func_str = "step(1.-sqrt(x))*(1.-sqrt(x))";
307 : } else {
308 : func_str = ktype;
309 : }
310 9 : std::string vm_str="";
311 9 : if( ktype=="von-misses" ) {
312 : vm_str=" VON_MISSES";
313 : }
314 :
315 9 : unsigned nvals = argnames.size();
316 : bool norm;
317 9 : parseFlag("NORMALIZED",norm);
318 9 : if( !usemahalanobis ) {
319 : // Invert the variance
320 4 : readInputLine( getShortcutLabel() + "_icov: CUSTOM ARG=" + getShortcutLabel() + "_var FUNC=1/x PERIODIC=NO");
321 : // Compute the distance between the center of the basin and the current configuration
322 10 : readInputLine( getShortcutLabel() + "_dist_2: NORMALIZED_EUCLIDEAN_DISTANCE SQUARED" + vm_str +" ARG1=" + argstr + " ARG2=" + refpoint + " METRIC=" + getShortcutLabel() + "_icov");
323 : // And compute a determinent for the input covariance matrix if it is required
324 2 : if( norm ) {
325 2 : if( ktype=="von-misses" ) {
326 0 : readInputLine( getShortcutLabel() + "_vec: CUSTOM ARG=" + getShortcutLabel() + "_icov FUNC=x PERIODIC=NO" );
327 : } else {
328 6 : readInputLine( getShortcutLabel() + "_det: PRODUCT ARG=" + getShortcutLabel() + "_var");
329 : }
330 : }
331 : } else {
332 : // Invert the input covariance matrix
333 14 : readInputLine( getShortcutLabel() + "_icov: INVERT_MATRIX ARG=" + getShortcutLabel() + "_cov" );
334 : // Compute the distance between the center of the basin and the current configuration
335 28 : readInputLine( getShortcutLabel() + "_dist_2: MAHALANOBIS_DISTANCE SQUARED ARG1=" + argstr + " ARG2=" + refpoint + " METRIC=" + getShortcutLabel() + "_icov " + vm_str );
336 : // And compute a determinent for the input covariance matrix if it is required
337 7 : if( norm ) {
338 7 : if( ktype=="von-misses" ) {
339 21 : readInputLine( getShortcutLabel() + "_det: DIAGONALIZE ARG=" + getShortcutLabel() + "_cov VECTORS=all" );
340 : std::string num;
341 7 : std::string thisargname= getShortcutLabel() + "_det.vals-1";
342 14 : for(unsigned i=1; i<nvals; ++i) {
343 7 : Tools::convert( i+1, num );
344 21 : thisargname += "," + getShortcutLabel() + "_det.vals-" + num;
345 : }
346 7 : readInputLine( getShortcutLabel() + "_comp: CONCATENATE ARG=" + thisargname );
347 21 : readInputLine( getShortcutLabel() + "_vec: CUSTOM ARG=" + getShortcutLabel() + "_comp FUNC=1/x PERIODIC=NO");
348 : } else {
349 0 : readInputLine( getShortcutLabel() + "_det: DETERMINANT ARG=" + getShortcutLabel() + "_cov");
350 : }
351 : }
352 : }
353 :
354 : // Compute the Gaussian
355 : std::string wstr;
356 9 : parse("WEIGHT",wstr);
357 9 : if( norm ) {
358 9 : if( ktype=="gaussian" ) {
359 : std::string pstr;
360 2 : Tools::convert( sqrt(pow(2*pi,nvals)), pstr );
361 8 : readInputLine( getShortcutLabel() + "_vol: CUSTOM ARG=" + getShortcutLabel() + "_det FUNC=(sqrt(x)*" + pstr + ") PERIODIC=NO");
362 7 : } else if( ktype=="von-misses" ) {
363 : std::string min;
364 : std::string max;
365 7 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_dist_2_diff" );
366 7 : plumed_assert( av );
367 7 : if( !av->copyOutput(0)->isPeriodic() ) {
368 0 : error("VON_MISSES only works with periodic variables");
369 : }
370 7 : av->copyOutput(0)->getDomain(min,max);
371 14 : readInputLine( getShortcutLabel() + "_bes: BESSEL ORDER=0 ARG=" + getShortcutLabel() + "_vec");
372 28 : readInputLine( getShortcutLabel() + "_cc: CUSTOM ARG=" + getShortcutLabel() + "_bes FUNC=("+max+"-"+min+")*x PERIODIC=NO");
373 21 : readInputLine( getShortcutLabel() + "_vol: PRODUCT ARG=" + getShortcutLabel() + "_cc");
374 : } else {
375 0 : error("only gaussian and von-misses kernels are normalizable");
376 : }
377 : // And the (suitably normalized) kernel
378 45 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_dist_2," + getShortcutLabel() + "_vol FUNC=" + wstr + "*exp(-x/2)/y PERIODIC=NO");
379 : } else {
380 0 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_dist_2 FUNC=" + wstr + "*" + func_str + " PERIODIC=NO");
381 : }
382 9 : checkRead();
383 :
384 93 : }
385 :
386 : }
387 : }
388 :
389 :
|
