Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #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 : namespace PLMD {
38 :
39 : template<class T>
40 17000148 : bool Tools::convertToAny(const std::string & str,T & t) {
41 31768385 : std::istringstream istr(str.c_str());
42 17000148 : bool ok=static_cast<bool>(istr>>t);
43 17000148 : if(!ok) {
44 : return false;
45 : }
46 : std::string remaining;
47 14896401 : istr>>remaining;
48 14896401 : return remaining.length()==0;
49 17000148 : }
50 :
51 2228374 : bool Tools::convertNoexcept(const std::string & str,int & t) {
52 2228374 : return convertToInt(str,t);
53 : }
54 :
55 3 : bool Tools::convertNoexcept(const std::string & str,long int & t) {
56 3 : return convertToInt(str,t);
57 : }
58 :
59 87 : bool Tools::convertNoexcept(const std::string & str,long long int & t) {
60 87 : return convertToInt(str,t);
61 : }
62 :
63 5256 : bool Tools::convertNoexcept(const std::string & str,unsigned & t) {
64 5256 : return convertToInt(str,t);
65 : }
66 :
67 97 : bool Tools::convertNoexcept(const std::string & str,long unsigned & t) {
68 97 : return convertToInt(str,t);
69 : }
70 :
71 54 : bool Tools::convertNoexcept(const std::string & str,long long unsigned & t) {
72 54 : return convertToInt(str,t);
73 : }
74 :
75 455798 : bool Tools::convertNoexcept(const std::string & str,AtomNumber &a) {
76 : // Note: AtomNumber's are NOT converted as int, so as to
77 : // avoid using lepton conversions.
78 : unsigned i;
79 455798 : bool r=convertToAny(str,i);
80 455798 : if(r) {
81 447412 : a.setSerial(i);
82 : }
83 455798 : return r;
84 : }
85 :
86 : template<class T>
87 2233871 : bool Tools::convertToInt(const std::string & str,T & t) {
88 : // First try standard conversion
89 2233871 : if(convertToAny(str,t)) {
90 : return true;
91 : }
92 : // Then use lepton
93 : try {
94 2008031 : double r=lepton::Parser::parse(str).evaluate(lepton::Constants());
95 :
96 : // now sanity checks on the resulting number
97 :
98 : // it should not overflow the requested int type:
99 : // (see https://stackoverflow.com/a/526092)
100 2008031 : if(r>std::nextafter(std::numeric_limits<T>::max(), 0)) {
101 : return false;
102 : }
103 2008031 : if(r<std::nextafter(std::numeric_limits<T>::min(), 0)) {
104 : return false;
105 : }
106 :
107 : // do the actual conversion
108 2008031 : auto tmp=static_cast<T>(std::round(r));
109 :
110 : // it should be *very close* to itself if converted back to double
111 2008031 : double diff= r-static_cast<double>(tmp);
112 2008031 : if(diff*diff > 1e-20) {
113 : return false;
114 : }
115 : // this is to accomodate small numerical errors and allow e.g. exp(log(7)) to be integer
116 :
117 : // it should be change if incremented or decremented by one (see https://stackoverflow.com/a/43656140)
118 2008027 : if(r == static_cast<double>(tmp-1)) {
119 : return false;
120 : }
121 2008026 : if(r == static_cast<double>(tmp+1)) {
122 : return false;
123 : }
124 :
125 : // everything is fine, then store in t
126 2008026 : t=tmp;
127 2008026 : return true;
128 0 : } catch(const PLMD::lepton::Exception& exc) {
129 : }
130 0 : return false;
131 : }
132 :
133 :
134 : template<class T>
135 14306942 : bool Tools::convertToReal(const std::string & str,T & t) {
136 14306942 : if(convertToAny(str,t)) {
137 : return true;
138 : }
139 8378270 : if(str=="PI" || str=="+PI" || str=="+pi" || str=="pi") {
140 1047323 : t=pi;
141 1047323 : return true;
142 2094614 : } else if(str=="-PI" || str=="-pi") {
143 1047242 : t=-pi;
144 1047242 : return true;
145 : }
146 : try {
147 65 : t=lepton::Parser::parse(str).evaluate(lepton::Constants());
148 24 : return true;
149 41 : } catch(const PLMD::lepton::Exception& exc) {
150 : }
151 41 : if( str.find("PI")!=std::string::npos ) {
152 0 : std::size_t pi_start=str.find_first_of("PI");
153 0 : if(str.substr(pi_start)!="PI") {
154 : return false;
155 : }
156 0 : std::istringstream nstr(str.substr(0,pi_start));
157 0 : T ff=0.0;
158 0 : bool ok=static_cast<bool>(nstr>>ff);
159 0 : if(!ok) {
160 : return false;
161 : }
162 0 : t=ff*pi;
163 : std::string remains;
164 0 : nstr>>remains;
165 0 : return remains.length()==0;
166 41 : } else if( str.find("pi")!=std::string::npos ) {
167 15 : std::size_t pi_start=str.find_first_of("pi");
168 30 : if(str.substr(pi_start)!="pi") {
169 : return false;
170 : }
171 15 : std::istringstream nstr(str.substr(0,pi_start));
172 15 : T ff=0.0;
173 15 : bool ok=static_cast<bool>(nstr>>ff);
174 15 : if(!ok) {
175 : return false;
176 : }
177 15 : t=ff*pi;
178 : std::string remains;
179 15 : nstr>>remains;
180 15 : return remains.length()==0;
181 41 : } else if(str=="NAN") {
182 0 : t=std::numeric_limits<double>::quiet_NaN();
183 0 : return true;
184 : }
185 : return false;
186 : }
187 :
188 0 : bool Tools::convertNoexcept(const std::string & str,float & t) {
189 0 : return convertToReal(str,t);
190 : }
191 :
192 14306814 : bool Tools::convertNoexcept(const std::string & str,double & t) {
193 14306814 : return convertToReal(str,t);
194 : }
195 :
196 128 : bool Tools::convertNoexcept(const std::string & str,long double & t) {
197 128 : return convertToReal(str,t);
198 : }
199 :
200 81337 : bool Tools::convertNoexcept(const std::string & str,std::string & t) {
201 : t=str;
202 81337 : return true;
203 : }
204 :
205 3912351 : std::vector<std::string> Tools::getWords(const std::string & line,const char* separators,int * parlevel,const char* parenthesis, const bool& delete_parenthesis) {
206 3912351 : plumed_massert(std::strlen(parenthesis)==1,"multiple parenthesis type not available");
207 3912351 : plumed_massert(parenthesis[0]=='(' || parenthesis[0]=='[' || parenthesis[0]=='{',
208 : "only ( [ { allowed as parenthesis");
209 3912351 : if(!separators) {
210 : separators=" \t\n";
211 : }
212 3912351 : const std::string sep(separators);
213 3912351 : char openpar=parenthesis[0];
214 : char closepar;
215 : if(openpar=='(') {
216 : closepar=')';
217 : }
218 3912351 : if(openpar=='[') {
219 : closepar=']';
220 : }
221 3912351 : if(openpar=='{') {
222 : closepar='}';
223 : }
224 : std::vector<std::string> words;
225 : std::string word;
226 : int parenthesisLevel=0;
227 3912351 : if(parlevel) {
228 26822 : parenthesisLevel=*parlevel;
229 : }
230 201059110 : for(unsigned i=0; i<line.length(); i++) {
231 : bool found=false;
232 : bool onParenthesis=false;
233 197146759 : if( (line[i]==openpar || line[i]==closepar) && delete_parenthesis ) {
234 : onParenthesis=true;
235 : }
236 197146759 : if(line[i]==closepar) {
237 4569 : parenthesisLevel--;
238 4569 : plumed_massert(parenthesisLevel>=0,"Extra closed parenthesis in '" + line + "'");
239 : }
240 197146759 : if(parenthesisLevel==0)
241 788666002 : for(unsigned j=0; j<sep.length(); j++)
242 591592385 : if(line[i]==sep[j]) {
243 : found=true;
244 : }
245 : // If at parenthesis level zero (outer)
246 197146759 : if(!(parenthesisLevel==0 && (found||onParenthesis))) {
247 140228247 : word.push_back(line[i]);
248 : }
249 : //if(onParenthesis) word.push_back(' ');
250 197146759 : if(line[i]==openpar) {
251 4569 : parenthesisLevel++;
252 : }
253 197146759 : if(found && word.length()>0) {
254 11818787 : if(!parlevel) {
255 11773604 : plumed_massert(parenthesisLevel==0,"Unmatching parenthesis in '" + line + "'");
256 : }
257 11818787 : words.push_back(word);
258 : word.clear();
259 : }
260 : }
261 3912351 : if(word.length()>0) {
262 3792872 : if(!parlevel) {
263 3766296 : plumed_massert(parenthesisLevel==0,"Unmatching parenthesis in '" + line + "'");
264 : }
265 3792872 : words.push_back(word);
266 : }
267 3912351 : if(parlevel) {
268 26822 : *parlevel=parenthesisLevel;
269 : }
270 3912351 : return words;
271 0 : }
272 :
273 15273 : bool Tools::getParsedLine(IFile& ifile,std::vector<std::string> & words, bool trimcomments) {
274 15273 : std::string line("");
275 : words.clear();
276 : bool stat;
277 : bool inside=false;
278 15273 : int parlevel=0;
279 : bool mergenext=false;
280 33850 : while((stat=ifile.getline(line))) {
281 32953 : if(trimcomments) {
282 32953 : trimComments(line);
283 : }
284 32953 : trim(line);
285 32953 : if(line.length()==0) {
286 6131 : continue;
287 : }
288 26822 : std::vector<std::string> w=getWords(line,NULL,&parlevel,"{",trimcomments);
289 26822 : if(!w.empty()) {
290 38932 : if(inside && *(w.begin())=="...") {
291 : inside=false;
292 1281 : if(w.size()==2) {
293 1126 : 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]);
294 : }
295 1281 : plumed_massert(w.size()<=2,"terminating \"...\" lines cannot consist of more than two words");
296 : w.clear();
297 1281 : if(!trimcomments) {
298 0 : words.push_back("...");
299 : }
300 25295 : } else if(*(w.end()-1)=="...") {
301 : inside=true;
302 : w.erase(w.end()-1);
303 : };
304 : int i0=0;
305 26576 : if(mergenext && words.size()>0 && w.size()>0) {
306 128 : words[words.size()-1]+=" "+w[0];
307 : i0=1;
308 : }
309 94583 : for(unsigned i=i0; i<w.size(); ++i) {
310 68007 : words.push_back(w[i]);
311 : }
312 : }
313 26822 : mergenext=(parlevel>0);
314 26822 : if(!inside) {
315 : break;
316 : }
317 12446 : if(!trimcomments && parlevel==0) {
318 0 : words.push_back("@newline");
319 12446 : } else if(!trimcomments) {
320 : words[words.size()-1] += " @newline";
321 : }
322 26822 : }
323 15273 : plumed_massert(parlevel==0,"non matching parenthesis");
324 15273 : if(words.size()>0) {
325 14220 : return true;
326 : }
327 : return stat;
328 : }
329 :
330 :
331 2804049 : bool Tools::getline(FILE* fp,std::string & line) {
332 : line="";
333 : const int bufferlength=1024;
334 : char buffer[bufferlength];
335 : bool ret;
336 2874150225 : for(int i=0; i<bufferlength; i++) {
337 2871346176 : buffer[i]='\0';
338 : }
339 2804049 : while((ret=fgets(buffer,bufferlength,fp))) {
340 2803302 : line.append(buffer);
341 2803302 : unsigned ss=std::strlen(buffer);
342 2803302 : if(ss>0)
343 2803302 : if(buffer[ss-1]=='\n') {
344 : break;
345 : }
346 : };
347 2804049 : if(line.length()>0)
348 2803302 : if(*(line.end()-1)=='\n') {
349 2803302 : line.erase(line.end()-1);
350 : }
351 2804049 : if(line.length()>0)
352 2803290 : if(*(line.end()-1)=='\r') {
353 1180 : line.erase(line.end()-1);
354 : }
355 2804049 : return ret;
356 : }
357 :
358 996441 : void Tools::trim(std::string & s) {
359 996441 : auto n=s.find_last_not_of(" \t");
360 996441 : if(n!=std::string::npos) {
361 983413 : s.resize(n+1);
362 : }
363 996441 : }
364 :
365 1356982 : void Tools::trimComments(std::string & s) {
366 1356982 : auto n=s.find_first_of("#");
367 1356982 : if(n!=std::string::npos) {
368 : s.resize(n);
369 : }
370 1356982 : }
371 :
372 383090 : bool Tools::caseInSensStringCompare(const std::string & str1, const std::string &str2) {
373 383090 : return ((str1.size() == str2.size()) && std::equal(str1.begin(), str1.end(), str2.begin(), [](char c1, char c2) {
374 712818 : return (c1 == c2 || std::toupper(c1) == std::toupper(c2));
375 383090 : }));
376 : }
377 :
378 112458 : bool Tools::getKey(std::vector<std::string>& line,const std::string & key,std::string & s,int rep) {
379 : s.clear();
380 437321 : for(auto p=line.begin(); p!=line.end(); ++p) {
381 383090 : if((*p).length()==0) {
382 0 : continue;
383 : }
384 383090 : std::string x=(*p).substr(0,key.length());
385 383090 : if(caseInSensStringCompare(x,key)) {
386 58227 : if((*p).length()==key.length()) {
387 : return false;
388 : }
389 58226 : std::string tmp=(*p).substr(key.length(),(*p).length());
390 : line.erase(p);
391 : s=tmp;
392 58226 : const std::string multi("@replicas:");
393 58226 : if(rep>=0 && startWith(s,multi)) {
394 24 : s=s.substr(multi.length(),s.length());
395 24 : std::vector<std::string> words=getWords(s,"\t\n ,");
396 24 : plumed_massert(rep<static_cast<int>(words.size()),"Number of fields in " + s + " not consistent with number of replicas");
397 24 : s=words[rep];
398 24 : }
399 : return true;
400 : }
401 : };
402 : return false;
403 : }
404 :
405 11965 : void Tools::interpretRanges(std::vector<std::string>&s) {
406 : std::vector<std::string> news;
407 58979 : for(const auto & p :s) {
408 47014 : news.push_back(p);
409 47014 : size_t dash=p.find("-");
410 47014 : if(dash==std::string::npos) {
411 45640 : continue;
412 : }
413 : int first;
414 4254 : if(!Tools::convertToAny(p.substr(0,dash),first)) {
415 753 : continue;
416 : }
417 1374 : int stride=1;
418 : int second;
419 1374 : size_t colon=p.substr(dash+1).find(":");
420 1374 : if(colon!=std::string::npos) {
421 108 : if(!Tools::convertToAny(p.substr(dash+1).substr(0,colon),second) ||
422 144 : !Tools::convertToAny(p.substr(dash+1).substr(colon+1),stride)) {
423 0 : continue;
424 : }
425 : } else {
426 2676 : if(!Tools::convertToAny(p.substr(dash+1),second)) {
427 0 : continue;
428 : }
429 : }
430 1374 : news.resize(news.size()-1);
431 1374 : if(first<=second) {
432 1373 : plumed_massert(stride>0,"interpreting ranges "+ p + ", stride should be positive");
433 408015 : for(int i=first; i<=second; i+=stride) {
434 : std::string ss;
435 406642 : convert(i,ss);
436 406642 : news.push_back(ss);
437 : }
438 : } else {
439 1 : plumed_massert(stride<0,"interpreting ranges "+ p + ", stride should be positive");
440 3 : for(int i=first; i>=second; i+=stride) {
441 : std::string ss;
442 2 : convert(i,ss);
443 2 : news.push_back(ss);
444 : }
445 : }
446 : }
447 11965 : s=news;
448 11965 : }
449 :
450 14552 : void Tools::interpretLabel(std::vector<std::string>&s) {
451 14552 : if(s.size()<2) {
452 333 : return;
453 : }
454 14219 : std::string s0=s[0];
455 14219 : unsigned l=s0.length();
456 14219 : if(l<1) {
457 : return;
458 : }
459 14219 : if(s0[l-1]==':') {
460 : s[0]=s[1];
461 21702 : s[1]="LABEL="+s0.substr(0,l-1);
462 : }
463 14219 : std::transform(s[0].begin(), s[0].end(), s[0].begin(), ::toupper);
464 : }
465 :
466 7968 : std::vector<std::string> Tools::ls(const std::string&d) {
467 : DIR*dir;
468 : std::vector<std::string> result;
469 7968 : if ((dir=opendir(d.c_str()))) {
470 : #if defined(__PLUMED_HAS_READDIR_R)
471 : struct dirent ent;
472 : #endif
473 : while(true) {
474 : struct dirent *res;
475 : #if defined(__PLUMED_HAS_READDIR_R)
476 : readdir_r(dir,&ent,&res);
477 : #else
478 132416 : res=readdir(dir);
479 : #endif
480 132416 : if(!res) {
481 : break;
482 : }
483 489824 : if(std::string(res->d_name)!="." && std::string(res->d_name)!="..") {
484 217024 : result.push_back(res->d_name);
485 : }
486 : }
487 7968 : closedir (dir);
488 : }
489 7968 : return result;
490 0 : }
491 :
492 4452 : void Tools::stripLeadingAndTrailingBlanks( std::string& str ) {
493 4452 : std::size_t first=str.find_first_not_of(' ');
494 4452 : std::size_t last=str.find_last_not_of(' ');
495 4452 : if( first<=last && first!=std::string::npos) {
496 8830 : str=str.substr(first,last+1);
497 : }
498 4452 : }
499 :
500 12015 : std::string Tools::extension(const std::string&s) {
501 12015 : size_t n=s.find_last_of(".");
502 : std::string ext;
503 12015 : if(n!=std::string::npos && n+1<s.length() && n+5>=s.length()) {
504 8746 : ext=s.substr(n+1);
505 8746 : if(ext.find("/")!=std::string::npos) {
506 : ext="";
507 : }
508 8746 : std::string base=s.substr(0,n);
509 8746 : if(base.length()==0) {
510 : ext="";
511 : }
512 8746 : if(base.length()>0 && base[base.length()-1]=='/') {
513 : ext="";
514 : }
515 : }
516 12015 : return ext;
517 : }
518 :
519 14 : double Tools::bessel0( const double& val ) {
520 14 : if (std::abs(val)<3.75) {
521 2 : double y = Tools::fastpow( val/3.75, 2 );
522 2 : return 1 + y*(3.5156229 +y*(3.0899424 + y*(1.2067492+y*(0.2659732+y*(0.0360768+y*0.0045813)))));
523 : }
524 12 : double ax=std::abs(val), y=3.75/ax, bx=std::exp(ax)/std::sqrt(ax);
525 12 : 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)))))));
526 12 : return ax*bx;
527 : }
528 :
529 994566 : bool Tools::startWith(const std::string & full,const std::string &start) {
530 994566 : return (full.substr(0,start.length())==start);
531 : }
532 :
533 103332 : bool Tools::findKeyword(const std::vector<std::string>&line,const std::string&key) {
534 103332 : const std::string search(key+"=");
535 659203 : for(const auto & p : line) {
536 605365 : if(startWith(p,search)) {
537 : return true;
538 : }
539 : }
540 : return false;
541 : }
542 :
543 521 : Tools::DirectoryChanger::DirectoryChanger(const char*path) {
544 521 : if(!path) {
545 : return;
546 : }
547 521 : if(std::strlen(path)==0) {
548 : return;
549 : }
550 : #ifdef __PLUMED_HAS_GETCWD
551 2 : char* ret=getcwd(cwd,buffersize);
552 2 : plumed_assert(ret)<<"Name of current directory too long, increase buffer size";
553 : #else
554 : plumed_error()<<"You are trying to use DirectoryChanger but your system does not support getcwd";
555 : #endif
556 : #ifdef __PLUMED_HAS_CHDIR
557 2 : int r=chdir(path);
558 3 : plumed_assert(r==0) <<"Cannot chdir to directory "<<path<<". The directory must exist!";
559 : #else
560 : plumed_error()<<"You are trying to use DirectoryChanger but your system does not support chdir";
561 : #endif
562 : }
563 :
564 521 : Tools::DirectoryChanger::~DirectoryChanger() {
565 : #ifdef __PLUMED_HAS_CHDIR
566 520 : if(std::strlen(cwd)==0) {
567 519 : return;
568 : }
569 1 : int ret=chdir(cwd);
570 : // we cannot put an assertion here (in a destructor) otherwise cppcheck complains
571 : // we thus just report the problem
572 1 : if(ret!=0) {
573 0 : std::fprintf(stderr,"+++ WARNING: cannot cd back to directory %s\n",cwd);
574 : }
575 : #endif
576 520 : }
577 :
578 21147 : std::unique_ptr<std::lock_guard<std::mutex>> Tools::molfile_lock() {
579 : static std::mutex mtx;
580 21147 : return Tools::make_unique<std::lock_guard<std::mutex>>(mtx);
581 : }
582 :
583 : /// Internal tool, I am keeping it private for now
584 : namespace {
585 :
586 : class process_one_exception {
587 : std::string & msg;
588 : bool first=true;
589 134 : void update() {
590 134 : if(!first) {
591 14 : msg+="\n\nThe above exception was the direct cause of the following exception:\n";
592 : }
593 134 : first=false;
594 134 : }
595 : public:
596 120 : process_one_exception(std::string & msg):
597 120 : msg(msg)
598 : {}
599 133 : void operator()(const std::exception & e) {
600 133 : update();
601 133 : msg+=e.what();
602 133 : }
603 0 : void operator()(const std::string & e) {
604 0 : update();
605 0 : msg+=e;
606 0 : }
607 1 : void operator()(const char* e) {
608 1 : update();
609 1 : msg+=e;
610 1 : }
611 : };
612 :
613 : template<class T>
614 134 : static void process_all_exceptions(T&& f) {
615 : try {
616 : // First throw the current exception
617 134 : throw;
618 148 : } catch(const std::nested_exception & e) {
619 : // If nested, we go recursive
620 : // notice that we apply function f only if exception is also a std::exception
621 : try {
622 14 : e.rethrow_nested();
623 28 : } catch(...) {
624 14 : process_all_exceptions(f);
625 : }
626 14 : auto d=dynamic_cast<const std::exception*>(&e);
627 14 : if(d) {
628 14 : f(*d);
629 : }
630 238 : } catch(const std::exception &e) {
631 : // If not nested, we end recursion
632 119 : f(e);
633 0 : } catch(const std::string &e) {
634 : // If not nested, we end recursion
635 0 : f(e);
636 2 : } catch(const char* e) {
637 : // If not nested, we end recursion
638 1 : f(e);
639 0 : } catch(...) {
640 : // If not nested and of unknown type, we stop the chain
641 : }
642 134 : }
643 :
644 : }
645 :
646 120 : std::string Tools::concatenateExceptionMessages() {
647 : std::string msg;
648 120 : process_all_exceptions(process_one_exception(msg));
649 120 : return msg;
650 : }
651 :
652 : }
|
