Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-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 "Random.h"
      23             : #include <cmath>
      24             : #include <cstdlib>
      25             : #include <sstream>
      26             : #include <iostream>
      27             : #include <iterator>
      28             : #include <functional>
      29             : 
      30             : namespace PLMD {
      31             : 
      32             : const double Random::fact=5.9604644775390625e-8;     /* 1 / 2^24  */
      33             : const double Random::EPS=3.0e-16;
      34             : const double Random::AM=1.0/IM;
      35             : const double Random::RNMX=(1.0-EPS); // 1.0-EPS;
      36             : const std::string Random::noname="noname";
      37             : 
      38      809094 : Random::Random(const std::string & name):
      39      809094 :   switchGaussian(false),
      40      809094 :   saveGaussian(0.0),
      41             : // this is required because if a Random object is created during
      42             : // initialization than Random::noname could still be initialized.
      43             : // In practice: without this it is not possible to declare
      44             : // a static Random object without enforcing the order of the
      45             : // static constructors.
      46      809094 :   name(&name!=&noname?name:"noname") {
      47      809094 :   iy=0;
      48    26700102 :   for(unsigned i=0; i<NTAB; i++) {
      49    25891008 :     iv[i]=0;
      50             :   }
      51      809094 :   setSeed(0);
      52      809094 : }
      53             : 
      54      810226 : void Random::setSeed(int idum_) {
      55      810226 :   if(idum_>0) {
      56         138 :     idum_=-idum_;
      57             :   }
      58      810226 :   idum=idum_;
      59      810226 :   incPrec=false;
      60      810226 : }
      61             : 
      62     1553429 : double Random::RandU01 () {
      63     1553429 :   if (incPrec) {
      64           1 :     return U01d();
      65             :   } else {
      66     1553428 :     return U01();
      67             :   }
      68             : }
      69             : 
      70             : 
      71           2 : double Random::U01d () {
      72             :   double u;
      73           2 :   u = U01();
      74           2 :   u += U01() * fact;
      75           2 :   return (u < 1.0) ? u : (u - 1.0);
      76             : }
      77             : 
      78     5570468 : double Random::U01() {
      79             :   int j,k;
      80             :   double temp;
      81     5570468 :   if (idum <= 0 || !iy) {
      82        1691 :     if (-idum < 1) {
      83         862 :       idum=1;
      84             :     } else {
      85         829 :       idum = -idum;
      86             :     }
      87       69331 :     for (j=NTAB+7; j>=0; j--) {
      88       67640 :       k=idum/IQ;
      89       67640 :       idum=IA*(idum-k*IQ)-IR*k;
      90       67640 :       if (idum < 0) {
      91         904 :         idum += IM;
      92             :       }
      93       67640 :       if (j < NTAB) {
      94       54112 :         iv[j] = idum;
      95             :       }
      96             :     }
      97        1691 :     iy=iv[0];
      98             :   }
      99     5570468 :   k=idum/IQ;
     100     5570468 :   idum=IA*(idum-k*IQ)-IR*k;
     101     5570468 :   if (idum < 0) {
     102       62498 :     idum += IM;
     103             :   }
     104     5570468 :   j=iy/NDIV;
     105     5570468 :   iy=iv[j];
     106     5570468 :   iv[j] = idum;
     107     5570468 :   if ((temp=AM*iy) > RNMX) {
     108             :     return RNMX;
     109             :   } else {
     110     5570468 :     return temp;
     111             :   }
     112             : }
     113             : 
     114           1 : void Random::WriteStateFull(std::ostream & out)const {
     115             :   out<<name<<std::endl;
     116           1 :   out<<idum<<" "<<iy;
     117          33 :   for(int i=0; i<NTAB; i++) {
     118          32 :     out<<" "<<iv[i];
     119             :   };
     120           1 :   out<<" "<<switchGaussian;
     121           1 :   out<<" "<<saveGaussian;
     122             :   out<<std::endl;
     123           1 : }
     124             : 
     125           0 : void Random::ReadStateFull (std::istream & in) {
     126           0 :   getline(in,name);
     127           0 :   in>>idum>>iy;
     128           0 :   for (int i = 0; i < NTAB; i++) {
     129           0 :     in>>iv[i];
     130             :   }
     131           0 :   in>>switchGaussian;
     132           0 :   in>>saveGaussian;
     133           0 : }
     134             : 
     135           2 : void Random::toString(std::string & str)const {
     136           2 :   std::ostringstream ostr;
     137           2 :   ostr<<idum<<"|"<<iy;
     138          66 :   for(int i=0; i<NTAB; i++) {
     139          64 :     ostr<<"|"<<iv[i];
     140             :   };
     141           2 :   str=ostr.str();
     142           2 : }
     143             : 
     144           2 : void Random::fromString(const std::string & str) {
     145           2 :   std::string s=str;
     146         425 :   for(unsigned i=0; i<s.length(); i++)
     147         423 :     if(s[i]=='|') {
     148          66 :       s[i]=' ';
     149             :     }
     150           2 :   std::istringstream istr(s.c_str());
     151           2 :   istr>>idum>>iy;
     152          66 :   for (int i = 0; i < NTAB; i++) {
     153          64 :     istr>>iv[i];
     154             :   }
     155           4 : }
     156             : 
     157             : // This allows to have the same stream of random numbers
     158             : // with different compilers:
     159             : #ifdef __INTEL_COMPILER
     160             : #pragma intel optimization_level 0
     161             : #endif
     162             : 
     163     1219531 : double Random::Gaussian() {
     164             :   double v1,v2,rsq;
     165     1219531 :   if(switchGaussian) {
     166      609739 :     switchGaussian=false;
     167      609739 :     return saveGaussian;
     168             :   }
     169             :   while(true) {
     170      776040 :     v1=2.0*RandU01()-1.0;
     171      776040 :     v2=2.0*RandU01()-1.0;
     172      776040 :     rsq=v1*v1+v2*v2;
     173      776040 :     if(rsq<1.0 && rsq>0.0) {
     174             :       break;
     175             :     }
     176             :   }
     177      609792 :   double fac=std::sqrt(-2.*std::log(rsq)/rsq);
     178      609792 :   saveGaussian=v1*fac;
     179      609792 :   switchGaussian=true;
     180      609792 :   return v2*fac;
     181             : }
     182             : 
     183           8 : void Random::Shuffle(std::vector<unsigned>& vec) {
     184             :   std::iterator_traits<std::vector<unsigned>::iterator >::difference_type i, n;
     185             :   n = vec.end() - vec.begin();
     186          12 :   for(i=n-1; i>0; --i) {
     187           4 :     std::swap(vec[i], vec[(int)round(RandU01() * IM) % i]);
     188             :   }
     189           8 : }
     190             : 
     191             : 
     192             : }