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 : #ifndef __PLUMED_tools_Tools_h
23 : #define __PLUMED_tools_Tools_h
24 :
25 : #include "AtomNumber.h"
26 : #include "Vector.h"
27 : #include "Tensor.h"
28 : #include <vector>
29 : #include <string>
30 : #include <cctype>
31 : #include <cstdio>
32 : #include <cmath>
33 : #include <limits>
34 : #include <algorithm>
35 : #include <sstream>
36 : #include <memory>
37 : #include <cstddef>
38 : #include <queue>
39 : #include <mutex>
40 :
41 : namespace PLMD {
42 :
43 : class IFile;
44 :
45 : /// \ingroup TOOLBOX
46 : /// Very small non-zero number
47 : const double epsilon(std::numeric_limits<double>::epsilon());
48 :
49 : /// \ingroup TOOLBOX
50 : /// Boltzman constant in kj/K
51 : const double kBoltzmann(0.0083144621);
52 :
53 : /// \ingroup TOOLBOX
54 : /// PI
55 : const double pi(3.141592653589793238462643383279502884197169399375105820974944592307);
56 :
57 : const double dp2cutoff(6.25);
58 :
59 : const double dp2cutoffA=1.00193418799744762399; // 1.0/(1-std::exp(-dp2cutoff));
60 : const double dp2cutoffB=-.00193418799744762399; // -std::exp(-dp2cutoff)/(1-std::exp(-dp2cutoff));
61 :
62 30397 : inline static bool dp2cutoffNoStretch() {
63 30397 : static const auto* res=std::getenv("PLUMED_DP2CUTOFF_NOSTRETCH");
64 30397 : return res;
65 : }
66 :
67 : /// \ingroup TOOLBOX
68 : /// Empty class which just contains several (static) tools
69 : class Tools {
70 : /// class to convert a string to a generic type T
71 : template<class T>
72 : static bool convertToAny(const std::string & str,T &t);
73 : /// class to convert a string to a real type T.
74 : /// T should be either float, double, or long double
75 : template<class T>
76 : static bool convertToReal(const std::string & str,T &t);
77 : /// class to convert a string to a int type T
78 : template<class T>
79 : static bool convertToInt(const std::string & str,T &t);
80 : public:
81 : /// Split the line in words using separators.
82 : /// It also take into account parenthesis. Outer parenthesis found are removed from
83 : /// output, and the text between them is considered as a single word. Only the
84 : /// outer parenthesis are processed, to allow nesting them.
85 : /// parlevel, if not NULL, is increased or decreased according to the number of opened/closed parenthesis
86 : static std::vector<std::string> getWords(const std::string & line,const char* sep=NULL,int* parlevel=NULL,const char* parenthesis="{", const bool& delete_parenthesis=true);
87 : /// Get a line from the file pointer ifile
88 : static bool getline(FILE*,std::string & line);
89 : /// Get a parsed line from the file pointer ifile
90 : /// This function already takes care of joining continued lines and splitting the
91 : /// resulting line into an array of words
92 : static bool getParsedLine(IFile&ifile,std::vector<std::string> & line, const bool trimcomments=true);
93 : /// compare two string in a case insensitive manner
94 : static bool caseInSensStringCompare(const std::string & str1, const std::string &str2);
95 : /// Convert a string to a double, reading it
96 : static bool convertNoexcept(const std::string & str,double & t);
97 : /// Convert a string to a long double, reading it
98 : static bool convertNoexcept(const std::string & str,long double & t);
99 : /// Convert a string to a float, reading it
100 : static bool convertNoexcept(const std::string & str,float & t);
101 : /// Convert a string to a int, reading it
102 : static bool convertNoexcept(const std::string & str,int & t);
103 : /// Convert a string to a long int, reading it
104 : static bool convertNoexcept(const std::string & str,long int & t);
105 : /// Convert a string to a long long int, reading it
106 : static bool convertNoexcept(const std::string & str,long long int & t);
107 : /// Convert a string to an unsigned int, reading it
108 : static bool convertNoexcept(const std::string & str,unsigned & t);
109 : /// Convert a string to a long unsigned int, reading it
110 : static bool convertNoexcept(const std::string & str,long unsigned & t);
111 : /// Convert a string to a long long unsigned int, reading it
112 : static bool convertNoexcept(const std::string & str,long long unsigned & t);
113 : /// Convert a string to a atom number, reading it
114 : static bool convertNoexcept(const std::string & str,AtomNumber & t);
115 : /// Convert a string to a string (i.e. copy)
116 : static bool convertNoexcept(const std::string & str,std::string & t);
117 : /// Convert anything into a string
118 : template<typename T>
119 : static bool convertNoexcept(T i,std::string & str);
120 : /// Convert anything into anything, throwing an exception in case there is an error
121 : /// Remove trailing blanks
122 : static void trim(std::string & s);
123 : /// Remove trailing comments
124 : static void trimComments(std::string & s);
125 : /// Apply pbc for a unitary cell
126 : static double pbc(double);
127 : /// Retrieve a key from a vector of options.
128 : /// It finds a key starting with "key=" or equal to "key" and copy the
129 : /// part after the = on s. E.g.:
130 : /// line.push_back("aa=xx");
131 : /// getKey(line,"aa",s);
132 : /// will set s="xx"
133 : static bool getKey(std::vector<std::string>& line,const std::string & key,std::string & s,int rep=-1);
134 : /// Find a keyword on the input line, eventually deleting it, and saving its value to val
135 : template <typename T,typename U>
136 14737602 : static void convert(const T & t,U & u) {
137 14738431 : plumed_assert(convertNoexcept(t,u)) <<"Error converting "<<t;
138 14737601 : }
139 : template <typename T>
140 : static bool parse(std::vector<std::string>&line,const std::string&key,T&val,int rep=-1);
141 : /// Find a keyword on the input line, eventually deleting it, and saving its value to a vector
142 : template <class T>
143 : static bool parseVector(std::vector<std::string>&line,const std::string&key,std::vector<T>&val,int rep=-1);
144 : /// Find a keyword without arguments on the input line
145 : static bool parseFlag(std::vector<std::string>&line,const std::string&key,bool&val);
146 : /// Find a keyword on the input line, just reporting if it exists or not
147 : static bool findKeyword(const std::vector<std::string>&line,const std::string&key);
148 : /// Interpret atom ranges
149 : static void interpretRanges(std::vector<std::string>&);
150 : /// Remove duplicates from a vector of type T
151 : template <typename T>
152 : static void removeDuplicates(std::vector<T>& vec);
153 : /// interpret ":" syntax for labels
154 : static void interpretLabel(std::vector<std::string>&s);
155 : /// list files in a directory
156 : static std::vector<std::string> ls(const std::string&);
157 : /// removes leading and trailing blanks from a string
158 : static void stripLeadingAndTrailingBlanks( std::string& str );
159 : /// Extract the extensions from a file name.
160 : /// E.g.: extension("pippo.xyz")="xyz".
161 : /// It only returns extensions with a length between 1 and 4
162 : /// E.g.: extension("pippo.12345")="" whereas extenion("pippo.1234")="1234";
163 : /// It is also smart enough to detect "/", so that
164 : /// extension("pippo/.t")="" whereas extension("pippo/a.t")="t"
165 : static std::string extension(const std::string&);
166 : /// Fast int power
167 : static double fastpow(double base,int exp);
168 : /// Modified 0th-order Bessel function of the first kind
169 : static double bessel0(const double& val);
170 : /// Check if a string full starts with string start.
171 : /// Same as full.find(start)==0
172 : static bool startWith(const std::string & full,const std::string &start);
173 : /**
174 : Tool to create a vector of raw pointers from a vector of unique_pointers (const version).
175 : Returning a vector is fast in C++11. It can be used in order to feed a vector<unique_ptr<T>>
176 : to a function that takes a vector<T*>.
177 : \verbatim
178 : // some function that takes a vec
179 : void func(std::vector<Data*> & vec);
180 : std::vector<std::unique_ptr<Data>> vec;
181 : // func(vec); // does not compile
182 : func(Tools::unique2raw(vec)); // compiles
183 : \endverbatim
184 : Notice that the conversion is fast but takes
185 : some time to allocate the new vector and copy the pointers. In case the function
186 : acting on the vector<T*> is very fast and we do not want to add significant overhead,
187 : it might be convenient to store a separate set of raw pointers.
188 : \verbatim
189 : // some function that takes a vec
190 : void func(std::vector<Data*> & vec);
191 : std::vector<std::unique_ptr<Data>> vec;
192 :
193 : // conversion done only once:
194 : auto vec_ptr=Tools::unique2raw(vec);
195 :
196 : for(int i=0;i<1000;i++){
197 : func(vec_ptr);
198 : }
199 : \endverbatim
200 : */
201 : template <typename T>
202 : static std::vector<T*> unique2raw(const std::vector<std::unique_ptr<T>>&);
203 : /// Tool to create a vector of raw pointers from a vector of unique_pointers.
204 : /// See the non const version.
205 : template <typename T>
206 : static std::vector<const T*> unique2raw(const std::vector<std::unique_ptr<const T>>&);
207 : /// Tiny class that changes directory and comes back when going out of scope.
208 : /// In case system calls to change dir are not available it throws an exception.
209 : /// \warning By construction, changing directory breaks thread safety! Use with care.
210 : class DirectoryChanger {
211 : static const std::size_t buffersize=4096;
212 : char cwd[buffersize]= {0};
213 : public:
214 : explicit DirectoryChanger(const char*path);
215 : ~DirectoryChanger();
216 : };
217 : /// Mimic C++14 std::make_unique
218 : template<class T> struct _Unique_if {
219 : typedef std::unique_ptr<T> _Single_object;
220 : };
221 : template<class T> struct _Unique_if<T[]> {
222 : typedef std::unique_ptr<T[]> _Unknown_bound;
223 : };
224 : template<class T, std::size_t N> struct _Unique_if<T[N]> {
225 : typedef void _Known_bound;
226 : };
227 : template<class T, class... Args>
228 : static typename _Unique_if<T>::_Single_object
229 5815131 : make_unique(Args&&... args) {
230 9430972 : return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
231 : }
232 : template<class T>
233 : static typename _Unique_if<T>::_Unknown_bound
234 189 : make_unique(std::size_t n) {
235 : typedef typename std::remove_extent<T>::type U;
236 9473 : return std::unique_ptr<T>(new U[n]());
237 : }
238 : template<class T, class... Args>
239 : static typename _Unique_if<T>::_Known_bound
240 : make_unique(Args&&...) = delete;
241 :
242 : static void set_to_zero(double*ptr,unsigned n) {
243 24392902 : for(unsigned i=0; i<n; i++) {
244 24190788 : ptr[i]=0.0;
245 : }
246 : }
247 :
248 : template<unsigned n>
249 226744 : static void set_to_zero(std::vector<VectorGeneric<n>> & vec) {
250 226744 : unsigned s=vec.size();
251 226744 : if(s==0) {
252 : return;
253 : }
254 202114 : set_to_zero(&vec[0][0],s*n);
255 : }
256 :
257 : template<unsigned n,unsigned m>
258 : static void set_to_zero(std::vector<TensorGeneric<n,m>> & vec) {
259 : unsigned s=vec.size();
260 : if(s==0) {
261 : return;
262 : }
263 : set_to_zero(&vec[0](0,0),s*n*m);
264 : }
265 :
266 :
267 :
268 :
269 : /// Merge sorted vectors.
270 : /// Takes a vector of pointers to containers and merge them.
271 : /// Containers should be already sorted.
272 : /// The content is appended to the result vector.
273 : /// Optionally, uses a priority_queue implementation.
274 : template<class C>
275 27378 : static void mergeSortedVectors(const std::vector<C*> & vecs, std::vector<typename C::value_type> & result,bool priority_queue=false) {
276 :
277 : /// local class storing the range of remaining objects to be pushed
278 : struct Entry {
279 : typename C::const_iterator fwdIt,endIt;
280 :
281 54587 : explicit Entry(C const& v) : fwdIt(v.begin()), endIt(v.end()) {}
282 : /// check if this vector still contains something to be pushed
283 : explicit operator bool () const {
284 : return fwdIt != endIt;
285 : }
286 : /// to allow using a priority_queu, which selects the highest element.
287 : /// we here (counterintuitively) define < as >
288 : bool operator< (Entry const& rhs) const {
289 : return *fwdIt > *rhs.fwdIt;
290 : }
291 : };
292 :
293 27378 : if(priority_queue) {
294 : std::priority_queue<Entry> queue;
295 : // note: queue does not have reserve() method
296 :
297 : // add vectors to the queue
298 : {
299 : std::size_t maxsize=0;
300 0 : for(unsigned i=0; i<vecs.size(); i++) {
301 0 : if(vecs[i]->size()>maxsize) {
302 : maxsize=vecs[i]->size();
303 : }
304 0 : if(!vecs[i]->empty()) {
305 0 : queue.push(Entry(*vecs[i]));
306 : }
307 : }
308 : // this is just to save multiple reallocations on push_back
309 0 : result.reserve(maxsize);
310 : }
311 :
312 : // first iteration (to avoid a if in the main loop)
313 0 : if(queue.empty()) {
314 : return;
315 : }
316 0 : auto tmp=queue.top();
317 0 : queue.pop();
318 0 : result.push_back(*tmp.fwdIt);
319 : tmp.fwdIt++;
320 0 : if(tmp) {
321 0 : queue.push(tmp);
322 : }
323 :
324 : // main loop
325 0 : while(!queue.empty()) {
326 0 : auto tmp=queue.top();
327 0 : queue.pop();
328 0 : if(result.back() < *tmp.fwdIt) {
329 0 : result.push_back(*tmp.fwdIt);
330 : }
331 : tmp.fwdIt++;
332 0 : if(tmp) {
333 0 : queue.push(tmp);
334 : }
335 : }
336 : } else {
337 :
338 : std::vector<Entry> entries;
339 27378 : entries.reserve(vecs.size());
340 :
341 : {
342 : std::size_t maxsize=0;
343 84706 : for(int i=0; i<vecs.size(); i++) {
344 57328 : if(vecs[i]->size()>maxsize) {
345 : maxsize=vecs[i]->size();
346 : }
347 57328 : if(!vecs[i]->empty()) {
348 54587 : entries.push_back(Entry(*vecs[i]));
349 : }
350 : }
351 : // this is just to save multiple reallocations on push_back
352 27378 : result.reserve(maxsize);
353 : }
354 :
355 424403 : while(!entries.empty()) {
356 : // find smallest pending element
357 : // we use max_element instead of min_element because we are defining < as > (see above)
358 397025 : const auto minval=*std::max_element(entries.begin(),entries.end())->fwdIt;
359 :
360 : // push it
361 397025 : result.push_back(minval);
362 :
363 : // fast forward vectors with elements equal to minval (to avoid duplicates)
364 1562323 : for(auto & e : entries)
365 1962677 : while(e && *e.fwdIt==minval) {
366 : ++e.fwdIt;
367 : }
368 :
369 : // remove from the entries vector all exhausted vectors
370 397025 : auto erase=std::remove_if(entries.begin(),entries.end(),[](const Entry & e) {
371 : return !e;
372 : });
373 : entries.erase(erase,entries.end());
374 : }
375 : }
376 :
377 : }
378 : static std::unique_ptr<std::lock_guard<std::mutex>> molfile_lock();
379 : /// Build a concatenated exception message.
380 : /// Should be called with an in-flight exception.
381 : static std::string concatenateExceptionMessages();
382 : };
383 :
384 : template <class T>
385 75820 : bool Tools::parse(std::vector<std::string>&line,const std::string&key,T&val,int rep) {
386 : std::string s;
387 151640 : if(!getKey(line,key+"=",s,rep)) {
388 : return false;
389 : }
390 31264 : if(s.length()>0 && !convertNoexcept(s,val)) {
391 3 : return false;
392 : }
393 : return true;
394 : }
395 :
396 : template <class T>
397 36605 : bool Tools::parseVector(std::vector<std::string>&line,const std::string&key,std::vector<T>&val,int rep) {
398 : std::string s;
399 73210 : if(!getKey(line,key+"=",s,rep)) {
400 : return false;
401 : }
402 : val.clear();
403 26929 : std::vector<std::string> words=getWords(s,"\t\n ,");
404 115814 : for(unsigned i=0; i<words.size(); ++i) {
405 : T v;
406 88885 : std::string s=words[i];
407 88885 : const std::string multi("@replicas:");
408 88885 : if(rep>=0 && startWith(s,multi)) {
409 6 : s=s.substr(multi.length(),s.length());
410 6 : std::vector<std::string> words=getWords(s,"\t\n ,");
411 6 : plumed_assert(rep<static_cast<int>(words.size()));
412 6 : s=words[rep];
413 6 : }
414 88885 : if(!convertNoexcept(s,v)) {
415 : return false;
416 : }
417 88885 : val.push_back(v);
418 : }
419 : return true;
420 26929 : }
421 :
422 : template<typename T>
423 13572 : void Tools::removeDuplicates(std::vector<T>& vec) {
424 13572 : std::sort(vec.begin(), vec.end());
425 13572 : vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
426 13572 : }
427 :
428 : inline
429 44485 : bool Tools::parseFlag(std::vector<std::string>&line,const std::string&key,bool&val) {
430 118481 : for(auto p=line.begin(); p!=line.end(); ++p) {
431 75628 : if(key==*p) {
432 1632 : val=true;
433 : line.erase(p);
434 : return true;
435 : }
436 : }
437 : return false;
438 : }
439 :
440 : /// beware: this brings any number into a pbc that ranges from -0.5 to 0.5
441 : inline
442 1856147440 : double Tools::pbc(double x) {
443 : #ifdef __PLUMED_PBC_WHILE
444 : while (x>0.5) {
445 : x-=1.0;
446 : }
447 : while (x<-0.5) {
448 : x+=1.0;
449 : }
450 : return x;
451 : #else
452 : if(std::numeric_limits<int>::round_style == std::round_toward_zero) {
453 : const double offset=100.0;
454 1856147440 : const double y=x+offset;
455 1856147440 : if(y>=0) {
456 1856142863 : return y-int(y+0.5);
457 : } else {
458 4577 : return y-int(y-0.5);
459 : }
460 : } else if(std::numeric_limits<int>::round_style == std::round_to_nearest) {
461 : return x-int(x);
462 : } else {
463 : return x-floor(x+0.5);
464 : }
465 : #endif
466 : }
467 :
468 : template<typename T>
469 540013 : bool Tools::convertNoexcept(T i,std::string & str) {
470 540013 : std::ostringstream ostr;
471 421822 : ostr<<i;
472 540013 : str=ostr.str();
473 540013 : return true;
474 540013 : }
475 :
476 : inline
477 79275828 : double Tools::fastpow(double base, int exp) {
478 79275828 : if(exp<0) {
479 0 : exp=-exp;
480 0 : base=1.0/base;
481 : }
482 : double result = 1.0;
483 375463658 : while (exp) {
484 262587636 : if (exp & 1) {
485 181434484 : result *= base;
486 : }
487 262587636 : exp >>= 1;
488 262587636 : base *= base;
489 : }
490 :
491 79275828 : return result;
492 : }
493 :
494 : template<typename T>
495 13703470 : std::vector<T*> Tools::unique2raw(const std::vector<std::unique_ptr<T>> & x) {
496 13703470 : std::vector<T*> v(x.size());
497 52721172 : for(unsigned i=0; i<x.size(); i++) {
498 39017702 : v[i]=x[i].get();
499 : }
500 13703470 : return v;
501 : }
502 :
503 : template<typename T>
504 : std::vector<const T*> Tools::unique2raw(const std::vector<std::unique_ptr<const T>> & x) {
505 : std::vector<const T*> v(x.size());
506 : for(unsigned i=0; i<x.size(); i++) {
507 : v[i]=x[i].get();
508 : }
509 : return v;
510 : }
511 :
512 : }
513 :
514 : #endif
515 :
|
