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             : 
      29             : namespace PLMD {
      30             : 
      31             : const std::string Random::noname="noname";
      32             : 
      33      813011 : Random::Random(const std::string & title):
      34      813011 :   switchGaussian(false),
      35      813011 :   saveGaussian(0.0),
      36             : // this is required because if a Random object is created during
      37             : // initialization than Random::noname could still be initialized.
      38             : // In practice: without this it is not possible to declare
      39             : // a static Random object without enforcing the order of the
      40             : // static constructors.
      41      813011 :   name(&title!=&noname?title:"noname") {
      42      813011 :   iy=0;
      43    26829363 :   for(unsigned i=0; i<NTAB; i++) {
      44    26016352 :     iv[i]=0;
      45             :   }
      46      813011 :   setSeed(0);
      47      813011 : }
      48             : 
      49      814609 : void Random::setSeed(int idum_) {
      50      814609 :   if(idum_>0) {
      51         362 :     idum_=-idum_;
      52             :   }
      53      814609 :   idum=idum_;
      54      814609 :   incPrec=false;
      55      814609 : }
      56             : 
      57     3751378 : double Random::RandU01 () {
      58     3751378 :   if (incPrec) {
      59           1 :     return U01d();
      60             :   } else {
      61     3751377 :     return U01();
      62             :   }
      63             : }
      64             : 
      65             : 
      66           2 : double Random::U01d () {
      67             :   double u;
      68           2 :   u = U01();
      69           2 :   u += U01() * fact;
      70           2 :   return (u < 1.0) ? u : (u - 1.0);
      71             : }
      72             : 
      73     7788021 : double Random::U01() {
      74             :   int j,k;
      75             :   double temp;
      76     7788021 :   if (idum <= 0 || !iy) {
      77        2123 :     if (-idum < 1) {
      78        1128 :       idum=1;
      79             :     } else {
      80         995 :       idum = -idum;
      81             :     }
      82       87043 :     for (j=NTAB+7; j>=0; j--) {
      83       84920 :       k=idum/IQ;
      84       84920 :       idum=IA*(idum-k*IQ)-IR*k;
      85       84920 :       if (idum < 0) {
      86        1238 :         idum += IM;
      87             :       }
      88       84920 :       if (j < NTAB) {
      89       67936 :         iv[j] = idum;
      90             :       }
      91             :     }
      92        2123 :     iy=iv[0];
      93             :   }
      94     7788021 :   k=idum/IQ;
      95     7788021 :   idum=IA*(idum-k*IQ)-IR*k;
      96     7788021 :   if (idum < 0) {
      97       87119 :     idum += IM;
      98             :   }
      99     7788021 :   j=iy/NDIV;
     100     7788021 :   iy=iv[j];
     101     7788021 :   iv[j] = idum;
     102     7788021 :   if ((temp=AM*iy) > RNMX) {
     103             :     return RNMX;
     104             :   } else {
     105     7788021 :     return temp;
     106             :   }
     107             : }
     108             : 
     109           1 : void Random::WriteStateFull(std::ostream & out)const {
     110             :   out<<name<<std::endl;
     111           1 :   out<<idum<<" "<<iy;
     112          33 :   for(int i=0; i<NTAB; i++) {
     113          32 :     out<<" "<<iv[i];
     114             :   };
     115           1 :   out<<" "<<switchGaussian;
     116           1 :   out<<" "<<saveGaussian;
     117             :   out<<std::endl;
     118           1 : }
     119             : 
     120           0 : void Random::ReadStateFull (std::istream & in) {
     121           0 :   getline(in,name);
     122           0 :   in>>idum>>iy;
     123           0 :   for (int i = 0; i < NTAB; i++) {
     124           0 :     in>>iv[i];
     125             :   }
     126           0 :   in>>switchGaussian;
     127           0 :   in>>saveGaussian;
     128           0 : }
     129             : 
     130           2 : void Random::toString(std::string & str)const {
     131           2 :   std::ostringstream ostr;
     132           2 :   ostr<<idum<<"|"<<iy;
     133          66 :   for(int i=0; i<NTAB; i++) {
     134          64 :     ostr<<"|"<<iv[i];
     135             :   };
     136           2 :   str=ostr.str();
     137           2 : }
     138             : 
     139           2 : void Random::fromString(const std::string & str) {
     140           2 :   std::string s=str;
     141         425 :   for(unsigned i=0; i<s.length(); i++)
     142         423 :     if(s[i]=='|') {
     143          66 :       s[i]=' ';
     144             :     }
     145           2 :   std::istringstream istr(s.c_str());
     146           2 :   istr>>idum>>iy;
     147          66 :   for (int i = 0; i < NTAB; i++) {
     148          64 :     istr>>iv[i];
     149             :   }
     150           4 : }
     151             : 
     152             : // This allows to have the same stream of random numbers
     153             : // with different compilers:
     154             : #ifdef __INTEL_COMPILER
     155             : #pragma intel optimization_level 0
     156             : #endif
     157             : 
     158     1386820 : double Random::Gaussian() {
     159             :   double v1,v2,rsq;
     160     1386820 :   if(switchGaussian) {
     161      693379 :     switchGaussian=false;
     162      693379 :     return saveGaussian;
     163             :   }
     164             :   while(true) {
     165      882469 :     v1=2.0*RandU01()-1.0;
     166      882469 :     v2=2.0*RandU01()-1.0;
     167      882469 :     rsq=v1*v1+v2*v2;
     168      882469 :     if(rsq<1.0 && rsq>0.0) {
     169             :       break;
     170             :     }
     171             :   }
     172      693441 :   double fac=std::sqrt(-2.*std::log(rsq)/rsq);
     173      693441 :   saveGaussian=v1*fac;
     174      693441 :   switchGaussian=true;
     175      693441 :   return v2*fac;
     176             : }
     177             : 
     178      360068 : void Random::Shuffle(std::vector<unsigned>& vec) {
     179             :   std::iterator_traits<std::vector<unsigned>::iterator >::difference_type i, n;
     180             :   n = vec.end() - vec.begin();
     181     1080300 :   for(i=n-1; i>0; --i) {
     182             :     // Fisher-Yates: j is drawn from [0,i]. Drawing it from [0,i-1] (Sattolo's algorithm)
     183             :     // would only produce cyclic permutations
     184      720232 :     std::swap(vec[i], vec[(int)round(RandU01() * IM) % (i+1)]);
     185             :   }
     186      360068 : }
     187             : 
     188             : 
     189             : }