1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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 "Tools.h"
23 #include "AtomNumber.h"
24 #include "Exception.h"
25 #include "IFile.h"
26 #include "lepton/Lepton.h"
27 #include <cstring>
28 #include <dirent.h>
29 #include <iostream>
30 #include <map>
31 #if defined(__PLUMED_HAS_CHDIR) || defined(__PLUMED_HAS_GETCWD)
32 #include <unistd.h>
33 #endif
34
35 #include <iomanip>
36
37 using namespace std;
38 namespace PLMD {
39
40 template<class T>
convertToAny(const string & str,T & t)41 bool Tools::convertToAny(const string & str,T & t) {
42 istringstream istr(str.c_str());
43 bool ok=static_cast<bool>(istr>>t);
44 if(!ok) return false;
45 string remaining;
46 istr>>remaining;
47 return remaining.length()==0;
48 }
49
convert(const string & str,int & t)50 bool Tools::convert(const string & str,int & t) {
51 return convertToInt(str,t);
52 }
53
convert(const string & str,long int & t)54 bool Tools::convert(const string & str,long int & t) {
55 return convertToInt(str,t);
56 }
57
convert(const string & str,unsigned & t)58 bool Tools::convert(const string & str,unsigned & t) {
59 return convertToInt(str,t);
60 }
61
convert(const string & str,long unsigned & t)62 bool Tools::convert(const string & str,long unsigned & t) {
63 return convertToInt(str,t);
64 }
65
convert(const string & str,AtomNumber & a)66 bool Tools::convert(const string & str,AtomNumber &a) {
67 // Note: AtomNumber's are NOT converted as int, so as to
68 // avoid using lepton conversions.
69 unsigned i;
70 bool r=convertToAny(str,i);
71 if(r) a.setSerial(i);
72 return r;
73 }
74
75 template<class T>
convertToInt(const string & str,T & t)76 bool Tools::convertToInt(const string & str,T & t) {
77 // First try standard conversion
78 if(convertToAny(str,t)) return true;
79 // Then use lepton
80 try {
81 double r=lepton::Parser::parse(str).evaluate(lepton::Constants());
82
83 // now sanity checks on the resulting number
84
85 // it should not overflow the requested int type:
86 // (see https://stackoverflow.com/a/526092)
87 if(r>std::nextafter(std::numeric_limits<T>::max(), 0)) return false;
88 if(r<std::nextafter(std::numeric_limits<T>::min(), 0)) return false;
89
90 // do the actual conversion
91 auto tmp=static_cast<T>(std::round(r));
92
93 // it should be *very close* to itself if converted back to double
94 double diff= r-static_cast<double>(tmp);
95 if(diff*diff > 1e-20) return false;
96 // this is to accomodate small numerical errors and allow e.g. exp(log(7)) to be integer
97
98 // it should be change if incremented or decremented by one (see https://stackoverflow.com/a/43656140)
99 if(r == static_cast<double>(tmp-1)) return false;
100 if(r == static_cast<double>(tmp+1)) return false;
101
102 // everything is fine, then store in t
103 t=tmp;
104 return true;
105 } catch(PLMD::lepton::Exception& exc) {
106 }
107 return false;
108 }
109
110
111 template<class T>
convertToReal(const string & str,T & t)112 bool Tools::convertToReal(const string & str,T & t) {
113 if(convertToAny(str,t)) return true;
114 if(str=="PI" || str=="+PI" || str=="+pi" || str=="pi") {
115 t=pi; return true;
116 } else if(str=="-PI" || str=="-pi") {
117 t=-pi; return true;
118 }
119 try {
120 t=lepton::Parser::parse(str).evaluate(lepton::Constants());
121 return true;
122 } catch(const PLMD::lepton::Exception& exc) {
123 }
124 if( str.find("PI")!=std::string::npos ) {
125 std::size_t pi_start=str.find_first_of("PI");
126 if(str.substr(pi_start)!="PI") return false;
127 istringstream nstr(str.substr(0,pi_start));
128 T ff=0.0; bool ok=static_cast<bool>(nstr>>ff);
129 if(!ok) return false;
130 t=ff*pi;
131 std::string remains; nstr>>remains;
132 return remains.length()==0;
133 } else if( str.find("pi")!=std::string::npos ) {
134 std::size_t pi_start=str.find_first_of("pi");
135 if(str.substr(pi_start)!="pi") return false;
136 istringstream nstr(str.substr(0,pi_start));
137 T ff=0.0; bool ok=static_cast<bool>(nstr>>ff);
138 if(!ok) return false;
139 t=ff*pi;
140 std::string remains; nstr>>remains;
141 return remains.length()==0;
142 } else if(str=="NAN") {
143 t=std::numeric_limits<double>::quiet_NaN();
144 return true;
145 }
146 return false;
147 }
148
convert(const string & str,float & t)149 bool Tools::convert(const string & str,float & t) {
150 return convertToReal(str,t);
151 }
152
convert(const string & str,double & t)153 bool Tools::convert(const string & str,double & t) {
154 return convertToReal(str,t);
155 }
156
convert(const string & str,long double & t)157 bool Tools::convert(const string & str,long double & t) {
158 return convertToReal(str,t);
159 }
160
convert(const string & str,string & t)161 bool Tools::convert(const string & str,string & t) {
162 t=str;
163 return true;
164 }
165
getWords(const string & line,const char * separators,int * parlevel,const char * parenthesis,const bool & delete_parenthesis)166 vector<string> Tools::getWords(const string & line,const char* separators,int * parlevel,const char* parenthesis, const bool& delete_parenthesis) {
167 plumed_massert(strlen(parenthesis)==1,"multiple parenthesis type not available");
168 plumed_massert(parenthesis[0]=='(' || parenthesis[0]=='[' || parenthesis[0]=='{',
169 "only ( [ { allowed as parenthesis");
170 if(!separators) separators=" \t\n";
171 const string sep(separators);
172 char openpar=parenthesis[0];
173 char closepar;
174 if(openpar=='(') closepar=')';
175 if(openpar=='[') closepar=']';
176 if(openpar=='{') closepar='}';
177 vector<string> words;
178 string word;
179 int parenthesisLevel=0;
180 if(parlevel) parenthesisLevel=*parlevel;
181 for(unsigned i=0; i<line.length(); i++) {
182 bool found=false;
183 bool onParenthesis=false;
184 if( (line[i]==openpar || line[i]==closepar) && delete_parenthesis ) onParenthesis=true;
185 if(line[i]==closepar) {
186 parenthesisLevel--;
187 plumed_massert(parenthesisLevel>=0,"Extra closed parenthesis in '" + line + "'");
188 }
189 if(parenthesisLevel==0) for(unsigned j=0; j<sep.length(); j++) if(line[i]==sep[j]) found=true;
190 // If at parenthesis level zero (outer)
191 if(!(parenthesisLevel==0 && (found||onParenthesis))) word.push_back(line[i]);
192 //if(onParenthesis) word.push_back(' ');
193 if(line[i]==openpar) parenthesisLevel++;
194 if(found && word.length()>0) {
195 if(!parlevel) plumed_massert(parenthesisLevel==0,"Unmatching parenthesis in '" + line + "'");
196 words.push_back(word);
197 word.clear();
198 }
199 }
200 if(word.length()>0) {
201 if(!parlevel) plumed_massert(parenthesisLevel==0,"Unmatching parenthesis in '" + line + "'");
202 words.push_back(word);
203 }
204 if(parlevel) *parlevel=parenthesisLevel;
205 return words;
206 }
207
getParsedLine(IFile & ifile,vector<string> & words,bool trimcomments)208 bool Tools::getParsedLine(IFile& ifile,vector<string> & words, bool trimcomments) {
209 string line("");
210 words.clear();
211 bool stat;
212 bool inside=false;
213 int parlevel=0;
214 bool mergenext=false;
215 while((stat=ifile.getline(line))) {
216 if(trimcomments) trimComments(line);
217 trim(line);
218 if(line.length()==0) continue;
219 vector<string> w=getWords(line,NULL,&parlevel,"{",trimcomments);
220 if(!w.empty()) {
221 if(inside && *(w.begin())=="...") {
222 inside=false;
223 if(w.size()==2) plumed_massert(w[1]==words[0],"second word in terminating \"...\" "+w[1]+" line, if present, should be equal to first word of directive: "+words[0]);
224 plumed_massert(w.size()<=2,"terminating \"...\" lines cannot consist of more than two words");
225 w.clear(); if(!trimcomments) words.push_back("...");
226 } else if(*(w.end()-1)=="...") {
227 inside=true;
228 w.erase(w.end()-1);
229 };
230 int i0=0;
231 if(mergenext && words.size()>0 && w.size()>0) {
232 words[words.size()-1]+=" "+w[0];
233 i0=1;
234 }
235 for(unsigned i=i0; i<w.size(); ++i) words.push_back(w[i]);
236 }
237 mergenext=(parlevel>0);
238 if(!inside)break;
239 if(!trimcomments && parlevel==0) words.push_back("@newline");
240 else if(!trimcomments) words[words.size()-1] += " @newline";
241 }
242 plumed_massert(parlevel==0,"non matching parenthesis");
243 if(words.size()>0) return true;
244 return stat;
245 }
246
247
getline(FILE * fp,string & line)248 bool Tools::getline(FILE* fp,string & line) {
249 line="";
250 const int bufferlength=1024;
251 char buffer[bufferlength];
252 bool ret;
253 for(int i=0; i<bufferlength; i++) buffer[i]='\0';
254 while((ret=fgets(buffer,bufferlength,fp))) {
255 line.append(buffer);
256 unsigned ss=strlen(buffer);
257 if(ss>0) if(buffer[ss-1]=='\n') break;
258 };
259 if(line.length()>0) if(*(line.end()-1)=='\n') line.erase(line.end()-1);
260 if(line.length()>0) if(*(line.end()-1)=='\r') line.erase(line.end()-1);
261 return ret;
262 }
263
trim(string & s)264 void Tools::trim(string & s) {
265 size_t n=s.find_last_not_of(" \t");
266 s=s.substr(0,n+1);
267 }
268
trimComments(string & s)269 void Tools::trimComments(string & s) {
270 size_t n=s.find_first_of("#");
271 s=s.substr(0,n);
272 }
273
caseInSensStringCompare(const std::string & str1,const std::string & str2)274 bool Tools::caseInSensStringCompare(const std::string & str1, const std::string &str2)
275 {
276 return ((str1.size() == str2.size()) && std::equal(str1.begin(), str1.end(), str2.begin(), [](char c1, char c2) {
277 return (c1 == c2 || std::toupper(c1) == std::toupper(c2));
278 }));
279 }
280
getKey(vector<string> & line,const string & key,string & s,int rep)281 bool Tools::getKey(vector<string>& line,const string & key,string & s,int rep) {
282 s.clear();
283 for(auto p=line.begin(); p!=line.end(); ++p) {
284 if((*p).length()==0) continue;
285 string x=(*p).substr(0,key.length());
286 if(caseInSensStringCompare(x,key)) {
287 if((*p).length()==key.length())return false;
288 string tmp=(*p).substr(key.length(),(*p).length());
289 line.erase(p);
290 s=tmp;
291 const std::string multi("@replicas:");
292 if(rep>=0 && startWith(s,multi)) {
293 s=s.substr(multi.length(),s.length());
294 std::vector<std::string> words=getWords(s,"\t\n ,");
295 plumed_massert(rep<static_cast<int>(words.size()),"Number of fields in " + s + " not consistent with number of replicas");
296 s=words[rep];
297 }
298 return true;
299 }
300 };
301 return false;
302 }
303
interpretRanges(std::vector<std::string> & s)304 void Tools::interpretRanges(std::vector<std::string>&s) {
305 vector<string> news;
306 for(const auto & p :s) {
307 news.push_back(p);
308 size_t dash=p.find("-");
309 if(dash==string::npos) continue;
310 int first;
311 if(!Tools::convertToAny(p.substr(0,dash),first)) continue;
312 int stride=1;
313 int second;
314 size_t colon=p.substr(dash+1).find(":");
315 if(colon!=string::npos) {
316 if(!Tools::convertToAny(p.substr(dash+1).substr(0,colon),second) ||
317 !Tools::convertToAny(p.substr(dash+1).substr(colon+1),stride)) continue;
318 } else {
319 if(!Tools::convertToAny(p.substr(dash+1),second)) continue;
320 }
321 news.resize(news.size()-1);
322 if(first<=second) {
323 plumed_massert(stride>0,"interpreting ranges "+ p + ", stride should be positive");
324 for(int i=first; i<=second; i+=stride) {
325 string ss;
326 convert(i,ss);
327 news.push_back(ss);
328 }
329 } else {
330 plumed_massert(stride<0,"interpreting ranges "+ p + ", stride should be positive");
331 for(int i=first; i>=second; i+=stride) {
332 string ss;
333 convert(i,ss);
334 news.push_back(ss);
335 }
336 }
337 }
338 s=news;
339 }
340
interpretLabel(vector<string> & s)341 void Tools::interpretLabel(vector<string>&s) {
342 if(s.size()<2)return;
343 string s0=s[0];
344 unsigned l=s0.length();
345 if(l<1) return;
346 if(s0[l-1]==':') {
347 s[0]=s[1];
348 s[1]="LABEL="+s0.substr(0,l-1);
349 }
350 std::transform(s[0].begin(), s[0].end(), s[0].begin(), ::toupper);
351 }
352
ls(const string & d)353 vector<string> Tools::ls(const string&d) {
354 DIR*dir;
355 vector<string> result;
356 if ((dir=opendir(d.c_str()))) {
357 #if defined(__PLUMED_HAS_READDIR_R)
358 struct dirent ent;
359 #endif
360 while(true) {
361 struct dirent *res;
362 #if defined(__PLUMED_HAS_READDIR_R)
363 readdir_r(dir,&ent,&res);
364 #else
365 res=readdir(dir);
366 #endif
367 if(!res) break;
368 if(string(res->d_name)!="." && string(res->d_name)!="..") result.push_back(res->d_name);
369 }
370 closedir (dir);
371 }
372 return result;
373 }
374
stripLeadingAndTrailingBlanks(std::string & str)375 void Tools::stripLeadingAndTrailingBlanks( std::string& str ) {
376 std::size_t first=str.find_first_not_of(' ');
377 std::size_t last=str.find_last_not_of(' ');
378 if( first<=last && first!=std::string::npos) str=str.substr(first,last+1);
379 }
380
extension(const std::string & s)381 std::string Tools::extension(const std::string&s) {
382 size_t n=s.find_last_of(".");
383 std::string ext;
384 if(n!=std::string::npos && n+1<s.length() && n+5>=s.length()) {
385 ext=s.substr(n+1);
386 if(ext.find("/")!=std::string::npos) ext="";
387 string base=s.substr(0,n);
388 if(base.length()==0) ext="";
389 if(base.length()>0 && base[base.length()-1]=='/') ext="";
390 }
391 return ext;
392 }
393
bessel0(const double & val)394 double Tools::bessel0( const double& val ) {
395 if (fabs(val)<3.75) {
396 double y = Tools::fastpow( val/3.75, 2 );
397 return 1 + y*(3.5156229 +y*(3.0899424 + y*(1.2067492+y*(0.2659732+y*(0.0360768+y*0.0045813)))));
398 }
399 double ax=fabs(val), y=3.75/ax, bx=std::exp(ax)/sqrt(ax);
400 ax=0.39894228+y*(0.01328592+y*(0.00225319+y*(-0.00157565+y*(0.00916281+y*(-0.02057706+y*(0.02635537+y*(-0.01647633+y*0.00392377)))))));
401 return ax*bx;
402 }
403
startWith(const std::string & full,const std::string & start)404 bool Tools::startWith(const std::string & full,const std::string &start) {
405 return (full.substr(0,start.length())==start);
406 }
407
findKeyword(const std::vector<std::string> & line,const std::string & key)408 bool Tools::findKeyword(const std::vector<std::string>&line,const std::string&key) {
409 const std::string search(key+"=");
410 for(const auto & p : line) {
411 if(startWith(p,search)) return true;
412 }
413 return false;
414 }
415
DirectoryChanger(const char * path)416 Tools::DirectoryChanger::DirectoryChanger(const char*path) {
417 if(!path) return;
418 if(std::strlen(path)==0) return;
419 #ifdef __PLUMED_HAS_GETCWD
420 char* ret=getcwd(cwd,buffersize);
421 plumed_assert(ret)<<"Name of current directory too long, increase buffer size";
422 #else
423 plumed_error()<<"You are trying to use DirectoryChanger but your system does not support getcwd";
424 #endif
425 #ifdef __PLUMED_HAS_CHDIR
426 int r=chdir(path);
427 plumed_assert(r==0) <<"Cannot chdir to directory "<<path<<". The directory must exist!";
428 #else
429 plumed_error()<<"You are trying to use DirectoryChanger but your system does not support chdir";
430 #endif
431 }
432
~DirectoryChanger()433 Tools::DirectoryChanger::~DirectoryChanger() {
434 #ifdef __PLUMED_HAS_CHDIR
435 if(strlen(cwd)==0) return;
436 int ret=chdir(cwd);
437 // we cannot put an assertion here (in a destructor) otherwise cppcheck complains
438 // we thus just report the problem
439 if(ret!=0) fprintf(stderr,"+++ WARNING: cannot cd back to directory %s\n",cwd);
440 #endif
441 }
442
443 }
444