Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "ActionWithArguments.h"
23 : #include "ActionWithValue.h"
24 : #include "tools/PDB.h"
25 : #include "PlumedMain.h"
26 : #include "ActionSet.h"
27 : #include <iostream>
28 : #ifdef __PLUMED_HAS_CREGEX
29 : #include <cstring>
30 : #include <regex.h>
31 : #endif
32 :
33 : using namespace std;
34 : namespace PLMD {
35 :
36 2431 : void ActionWithArguments::registerKeywords(Keywords& keys) {
37 : keys.reserve("numbered","ARG","the input for this action is the scalar output from one or more other actions. The particular scalars that you will use "
38 : "are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates "
39 : "a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the "
40 : "scalars calculated by all the proceeding actions in the input file are taken. Some actions have multi-component outputs and "
41 : "each component of the output has a specific label. For example a \\ref DISTANCE action labelled dist may have three components "
42 : "x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*."
43 : "More information on the referencing of Actions can be found in the section of the manual on the PLUMED \\ref Syntax. "
44 : "Scalar values can also be "
45 : "referenced using POSIX regular expressions as detailed in the section on \\ref Regex. To use this feature you you must compile "
46 9724 : "PLUMED with the appropriate flag.");
47 2431 : }
48 :
49 2344 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
50 4688 : std::string def; vector<string> c; arg.clear(); parseVector(key,c);
51 2689 : if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ) {
52 0 : if( keywords.getDefaultValue(key,def) ) c.push_back( def );
53 2341 : else return;
54 : }
55 2344 : interpretArgumentList(c,arg);
56 : }
57 :
58 6 : bool ActionWithArguments::parseArgumentList(const std::string&key,int i,std::vector<Value*>&arg) {
59 : vector<string> c;
60 : arg.clear();
61 6 : if(parseNumberedVector(key,i,c)) {
62 6 : interpretArgumentList(c,arg);
63 : return true;
64 6 : } else return false;
65 : }
66 :
67 2446 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, std::vector<Value*>&arg) {
68 14298 : for(unsigned i=0; i<c.size(); i++) {
69 : // is a regex? then just interpret it. The signal is ()
70 4706 : if(!c[i].compare(0,1,"(")) {
71 206 : unsigned l=c[i].length();
72 206 : if(!c[i].compare(l-1,1,")")) {
73 : // start regex parsing
74 : #ifdef __PLUMED_HAS_CREGEX
75 : // take the string enclosed in quotes and put in round brackets
76 205 : std::string myregex=c[i];
77 205 : log<<" Evaluating regexp for this action: "<<myregex<<"\n";
78 :
79 : regex_t reg; // regular expression
80 :
81 205 : int errcode=regcomp(®, myregex.c_str(),REG_EXTENDED|REG_NEWLINE); // compile the regular expression
82 205 : if(errcode) {
83 : // one can check the errors asking to regerror
84 1 : size_t errbuf_size = regerror(errcode, ®, NULL, 0);
85 2 : std::vector<char> errbuf(errbuf_size);
86 1 : regerror(errcode, ®, errbuf.data(), errbuf_size);
87 4 : plumed_error()<<"Error parsing regular expression: "<<errbuf.data();
88 : }
89 :
90 : // call regfree when reg goes out of scope
91 204 : auto deleter=[](regex_t* r) { regfree(r); };
92 : std::unique_ptr<regex_t,decltype(deleter)> reg_deleter(®,deleter);
93 :
94 204 : plumed_massert(reg.re_nsub==1,"I can parse with only one subexpression");
95 : regmatch_t match;
96 : // select all the actions that have a value
97 408 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
98 204 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
99 2184 : for(unsigned j=0; j<all.size(); j++) {
100 990 : std::vector<std::string> ss=all[j]->getComponentsVector();
101 440266 : for(unsigned k=0; k<ss.size(); ++k) {
102 219638 : unsigned ll=strlen(ss[k].c_str())+1;
103 219638 : std::vector<char> str(ll);
104 : strcpy(&str[0],ss[k].c_str());
105 : const char *ppstr=&str[0];
106 219638 : if(!regexec(®, ppstr, reg.re_nsub, &match, 0)) {
107 : //log.printf(" Something matched with \"%s\" : ",ss[k].c_str());
108 27768 : do {
109 27768 : if (match.rm_so != -1) { /* The regex is matching part of a string */
110 27768 : size_t matchlen = match.rm_eo - match.rm_so;
111 27768 : std::vector<char> submatch(matchlen+1);
112 27768 : strncpy(submatch.data(), ppstr+match.rm_so, matchlen+1);
113 27768 : submatch[matchlen]='\0';
114 : //log.printf(" subpattern %s\n", submatch.data());
115 : // this is the match: try to see if it is a valid action
116 27768 : std::string putativeVal(submatch.data());
117 27768 : if( all[j]->exists(putativeVal) ) {
118 42562 : arg.push_back(all[j]->copyOutput(putativeVal));
119 : //log.printf(" Action %s added! \n",putativeVal.c_str());
120 : }
121 : };
122 27768 : ppstr += match.rm_eo; /* Restart from last match */
123 27768 : } while(!regexec(®,ppstr,reg.re_nsub,&match,0));
124 : }
125 : }
126 990 : }
127 : #else
128 : plumed_merror("Regexp support not compiled!");
129 : #endif
130 : } else {
131 3 : plumed_merror("did you want to use regexp to input arguments? enclose it between two round braces (...) with no spaces!");
132 : }
133 : } else {
134 : std::size_t dot=c[i].find_first_of('.');
135 4500 : string a=c[i].substr(0,dot);
136 9000 : string name=c[i].substr(dot+1);
137 4500 : if(c[i].find(".")!=string::npos) { // if it contains a dot:
138 1607 : if(a=="*" && name=="*") {
139 : // Take all values from all actions
140 2 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
141 1 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
142 17 : for(unsigned j=0; j<all.size(); j++) {
143 56 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
144 : }
145 1598 : } else if ( name=="*") {
146 : // Take all the values from an action with a specific name
147 950 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a);
148 475 : if(!action) {
149 1 : std::string str=" (hint! the actions with value in this ActionSet are: ";
150 0 : str+=plumed.getActionSet().getLabelList<ActionWithValue*>()+")";
151 0 : error("cannot find action named " + a + str);
152 : }
153 475 : if( action->getNumberOfComponents()==0 ) error("found " + a +".* indicating use all components calculated by action with label " + a + " but this action has no components");
154 14915 : for(int k=0; k<action->getNumberOfComponents(); ++k) arg.push_back(action->copyOutput(k));
155 1123 : } else if ( a=="*" ) {
156 : // Take components from all actions with a specific name
157 14 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
158 7 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
159 : unsigned nval=0;
160 55 : for(unsigned j=0; j<all.size(); j++) {
161 72 : std::string flab; flab=all[j]->getLabel() + "." + name;
162 42 : if( all[j]->exists(flab) ) { arg.push_back(all[j]->copyOutput(flab)); nval++; }
163 : }
164 7 : if(nval==0) error("found no actions with a component called " + name );
165 : } else {
166 : // Take values with a specific name
167 2232 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a);
168 1116 : if(!action) {
169 1 : std::string str=" (hint! the actions with value in this ActionSet are: ";
170 0 : str+=plumed.getActionSet().getLabelList<ActionWithValue*>()+")";
171 0 : error("cannot find action named " + a +str);
172 : }
173 1116 : if( !(action->exists(c[i])) ) {
174 1 : std::string str=" (hint! the components in this actions are: ";
175 0 : str+=action->getComponentsList()+")";
176 0 : error("action " + a + " has no component named " + name + str);
177 : } ;
178 2233 : arg.push_back(action->copyOutput(c[i]));
179 : }
180 : } else { // if it doesn't contain a dot
181 2901 : if(c[i]=="*") {
182 : // Take all values from all actions
183 136 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
184 68 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
185 888 : for(unsigned j=0; j<all.size(); j++) {
186 2716 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
187 : }
188 : } else {
189 5666 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(c[i]);
190 2833 : if(!action) {
191 2 : std::string str=" (hint! the actions with value in this ActionSet are: ";
192 4 : str+=plumed.getActionSet().getLabelList<ActionWithValue*>()+")";
193 4 : error("cannot find action named " + c[i] + str );
194 : }
195 2832 : if( !(action->exists(c[i])) ) {
196 1 : std::string str=" (hint! the components in this actions are: ";
197 0 : str+=action->getComponentsList()+")";
198 0 : error("action " + c[i] + " has no component named " + c[i] +str);
199 : };
200 5665 : arg.push_back(action->copyOutput(c[i]));
201 : }
202 : }
203 : }
204 : }
205 2443 : }
206 :
207 416 : void ActionWithArguments::expandArgKeywordInPDB( PDB& pdb ) {
208 416 : std::vector<std::string> arg_names = pdb.getArgumentNames();
209 416 : if( arg_names.size()>0 ) {
210 : std::vector<Value*> arg_vals;
211 80 : interpretArgumentList( arg_names, arg_vals );
212 416 : }
213 416 : }
214 :
215 3462 : void ActionWithArguments::requestArguments(const vector<Value*> &arg) {
216 3462 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
217 3462 : arguments=arg;
218 3462 : clearDependencies();
219 : std::string fullname,name;
220 57124 : for(unsigned i=0; i<arguments.size(); i++) {
221 25100 : fullname=arguments[i]->getName();
222 25100 : if(fullname.find(".")!=string::npos) {
223 : std::size_t dot=fullname.find_first_of('.');
224 37862 : name=fullname.substr(0,dot);
225 : } else {
226 : name=fullname;
227 : }
228 50200 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
229 25100 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
230 25100 : addDependency(action);
231 : }
232 3462 : }
233 :
234 2347 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
235 : Action(ao),
236 4694 : lockRequestArguments(false)
237 : {
238 4697 : if( keywords.exists("ARG") ) {
239 : vector<Value*> arg;
240 4518 : parseArgumentList("ARG",arg);
241 :
242 2256 : if(!arg.empty()) {
243 2142 : log.printf(" with arguments");
244 68415 : for(unsigned i=0; i<arg.size(); i++) log.printf(" %s",arg[i]->getName().c_str());
245 2142 : log.printf("\n");
246 : }
247 2256 : requestArguments(arg);
248 : }
249 2344 : }
250 :
251 58 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
252 58 : if(!a) {
253 58 : a=dynamic_cast<ActionWithValue*>(this);
254 58 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
255 : }
256 :
257 : const int nval=a->getNumberOfComponents();
258 58 : const int npar=arguments.size();
259 58 : std::vector<double> value (nval*npar);
260 161 : for(int i=0; i<npar; i++) {
261 206 : double arg0=arguments[i]->get();
262 103 : arguments[i]->set(arg0+sqrt(epsilon));
263 103 : a->calculate();
264 103 : arguments[i]->set(arg0);
265 1216 : for(int j=0; j<nval; j++) {
266 2432 : value[i*nval+j]=a->getOutputQuantity(j);
267 : }
268 : }
269 58 : a->calculate();
270 58 : a->clearDerivatives();
271 1086 : for(int j=0; j<nval; j++) {
272 1086 : Value* v=a->copyOutput(j);
273 1886 : if( v->hasDerivatives() ) for(int i=0; i<npar; i++) v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/sqrt(epsilon));
274 : }
275 58 : }
276 :
277 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
278 522 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
279 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
280 261 : const Value* v1=arguments[i];
281 261 : const Value* v2=arguments[j];
282 261 : return Value::projection(*v1,*v2);
283 : }
284 :
285 350 : void ActionWithArguments::addForcesOnArguments( const std::vector<double>& forces ) {
286 1196 : for(unsigned i=0; i<arguments.size(); ++i) arguments[i]->addForce( forces[i] );
287 350 : }
288 :
289 5874 : }
|
