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      406024 : Random::Random(const std::string & name):
      39      406024 :   switchGaussian(false),
      40      406024 :   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      406024 :   name(&name!=&noname?name:"noname")
      47             : {
      48      406024 :   iy=0;
      49    13398792 :   for(unsigned i=0; i<NTAB; i++) iv[i]=0;
      50      406024 :   setSeed(0);
      51      406024 : }
      52             : 
      53      406835 : void Random::setSeed(int idum_) {
      54      406835 :   if(idum_>0) idum_=-idum_;
      55      406835 :   idum=idum_;
      56      406835 :   incPrec=false;
      57      406835 : }
      58             : 
      59     1582371 : double Random::RandU01 ()
      60             : {
      61     1582371 :   if (incPrec)
      62           1 :     return U01d();
      63             :   else
      64     1582370 :     return U01();
      65             : }
      66             : 
      67             : 
      68           2 : double Random::U01d ()
      69             : {
      70             :   double u;
      71           2 :   u = U01();
      72           2 :   u += U01() * fact;
      73           2 :   return (u < 1.0) ? u : (u - 1.0);
      74             : }
      75             : 
      76     5575120 : double Random::U01() {
      77             :   int j,k;
      78             :   double temp;
      79     5575120 :   if (idum <= 0 || !iy) {
      80        1285 :     if (-idum < 1) idum=1;
      81         626 :     else idum = -idum;
      82       52685 :     for (j=NTAB+7; j>=0; j--) {
      83       51400 :       k=idum/IQ;
      84       51400 :       idum=IA*(idum-k*IQ)-IR*k;
      85       51400 :       if (idum < 0) idum += IM;
      86       51400 :       if (j < NTAB) iv[j] = idum;
      87             :     }
      88        1285 :     iy=iv[0];
      89             :   }
      90     5575120 :   k=idum/IQ;
      91     5575120 :   idum=IA*(idum-k*IQ)-IR*k;
      92     5575120 :   if (idum < 0) idum += IM;
      93     5575120 :   j=iy/NDIV;
      94     5575120 :   iy=iv[j];
      95     5575120 :   iv[j] = idum;
      96     5575120 :   if ((temp=AM*iy) > RNMX) return RNMX;
      97     5575120 :   else return temp;
      98             : }
      99             : 
     100           1 : void Random::WriteStateFull(std::ostream & out)const {
     101             :   out<<name<<std::endl;
     102           1 :   out<<idum<<" "<<iy;
     103          33 :   for(int i=0; i<NTAB; i++) {
     104          32 :     out<<" "<<iv[i];
     105             :   };
     106           1 :   out<<" "<<switchGaussian;
     107           1 :   out<<" "<<saveGaussian;
     108             :   out<<std::endl;
     109           1 : }
     110             : 
     111           0 : void Random::ReadStateFull (std::istream & in) {
     112           0 :   getline(in,name);
     113           0 :   in>>idum>>iy;
     114           0 :   for (int i = 0; i < NTAB; i++) in>>iv[i];
     115           0 :   in>>switchGaussian;
     116           0 :   in>>saveGaussian;
     117           0 : }
     118             : 
     119           2 : void Random::toString(std::string & str)const {
     120           2 :   std::ostringstream ostr;
     121           2 :   ostr<<idum<<"|"<<iy;
     122          66 :   for(int i=0; i<NTAB; i++) {
     123          64 :     ostr<<"|"<<iv[i];
     124             :   };
     125           2 :   str=ostr.str();
     126           2 : }
     127             : 
     128           2 : void Random::fromString(const std::string & str) {
     129           2 :   std::string s=str;
     130         425 :   for(unsigned i=0; i<s.length(); i++) if(s[i]=='|') s[i]=' ';
     131           2 :   std::istringstream istr(s.c_str());
     132           2 :   istr>>idum>>iy;
     133          66 :   for (int i = 0; i < NTAB; i++) istr>>iv[i];
     134           4 : }
     135             : 
     136             : // This allows to have the same stream of random numbers
     137             : // with different compilers:
     138             : #ifdef __INTEL_COMPILER
     139             : #pragma intel optimization_level 0
     140             : #endif
     141             : 
     142     1242377 : double Random::Gaussian() {
     143             :   double v1,v2,rsq;
     144     1242377 :   if(switchGaussian) {
     145      621162 :     switchGaussian=false;
     146      621162 :     return saveGaussian;
     147             :   }
     148             :   while(true) {
     149      790679 :     v1=2.0*RandU01()-1.0;
     150      790679 :     v2=2.0*RandU01()-1.0;
     151      790679 :     rsq=v1*v1+v2*v2;
     152      790679 :     if(rsq<1.0 && rsq>0.0) break;
     153             :   }
     154      621215 :   double fac=std::sqrt(-2.*std::log(rsq)/rsq);
     155      621215 :   saveGaussian=v1*fac;
     156      621215 :   switchGaussian=true;
     157      621215 :   return v2*fac;
     158             : }
     159             : 
     160           2 : void Random::Shuffle(std::vector<unsigned>& vec) {
     161             :   std::iterator_traits<std::vector<unsigned>::iterator >::difference_type i, n;
     162             :   n = vec.end() - vec.begin();
     163           6 :   for(i=n-1; i>0; --i) {
     164           4 :     std::swap(vec[i], vec[(int)round(RandU01() * IM) % i]);
     165             :   }
     166           2 : }
     167             : 
     168             : 
     169             : }