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