NAMD
Public Member Functions | List of all members
Random Class Reference

#include <Random.h>

Public Member Functions

 Random (void)
 
 Random (unsigned long seed)
 
void init (unsigned long seed)
 
void skip (void)
 
void split (int iStream, int numStreams)
 
BigReal uniform (void)
 
BigReal gaussian (void)
 
BigReal sum_of_squared_gaussians (int64_t num_gaussians)
 Return the sum of i.i.d. squared Gaussians. More...
 
Vector gaussian_vector (void)
 
long integer (void)
 
template<class Elem >
void reorder (Elem *a, int n)
 
void uniform_array_f (float *a, int n)
 
void gaussian_array_f (float *a, int n)
 

Detailed Description

Definition at line 37 of file Random.h.

Constructor & Destructor Documentation

Random::Random ( void  )
inline

Definition at line 50 of file Random.h.

References init(), and RAND48_SEED.

50  {
51  init(0);
52  rand48_seed = RAND48_SEED;
53  }
void init(unsigned long seed)
Definition: Random.h:61
#define RAND48_SEED
Definition: Random.h:32
Random::Random ( unsigned long  seed)
inline

Definition at line 56 of file Random.h.

References init().

56  {
57  init(seed);
58  }
void init(unsigned long seed)
Definition: Random.h:61

Member Function Documentation

BigReal Random::gaussian ( void  )
inline

Definition at line 116 of file Random.h.

References uniform().

Referenced by Controller::adaptTempUpdate(), gaussian_array_f(), gaussian_vector(), Controller::langevinPiston1(), Controller::langevinPiston2(), colvarproxy_namd::rand_gaussian(), Controller::stochRescaleCoefficient(), and sum_of_squared_gaussians().

116  {
117  BigReal fac, r, v1, v2;
118 
119  if (second_gaussian_waiting) {
120  second_gaussian_waiting = 0;
121  return second_gaussian;
122  } else {
123  r = 2.; // r >= 1.523e-8 ensures abs result < 6
124  while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
125  v1 = 2.0 * uniform() - 1.0;
126  v2 = 2.0 * uniform() - 1.0;
127  r = v1*v1 + v2*v2;
128  }
129  fac = sqrt(-2.0 * log(r)/r);
130  // now make the Box-Muller transformation to get two normally
131  // distributed random numbers. Save one and return the other.
132  second_gaussian_waiting = 1;
133  second_gaussian = v1 * fac;
134  return v2 * fac;
135  }
136  }
BigReal uniform(void)
Definition: Random.h:109
double BigReal
Definition: common.h:112
void Random::gaussian_array_f ( float *  a,
int  n 
)
inline

Fill length n array from standard Gaussian distribution.

Definition at line 244 of file Random.h.

References gaussian().

244  {
245  // for now just call gaussian()
246  // ultimately we will provide faster implementation
247  for (int i=0; i < n; i++) {
248  a[i] = (float)gaussian();
249  }
250  }
BigReal gaussian(void)
Definition: Random.h:116
Vector Random::gaussian_vector ( void  )
inline
void Random::init ( unsigned long  seed)
inline

Definition at line 61 of file Random.h.

References RAND48_ADD, RAND48_MULT, RAND48_SEED, and UINT64_LITERAL.

Referenced by Random().

61  {
62  second_gaussian = 0;
63  second_gaussian_waiting = 0;
64  rand48_seed = seed & UINT64_LITERAL(0x00000000ffffffff);
65  rand48_seed = rand48_seed << 16;
66  rand48_seed |= RAND48_SEED & UINT64_LITERAL(0x0000ffff);
67  rand48_mult = RAND48_MULT;
68  rand48_add = RAND48_ADD;
69  }
#define RAND48_ADD
Definition: Random.h:34
#define RAND48_SEED
Definition: Random.h:32
#define RAND48_MULT
Definition: Random.h:33
#define UINT64_LITERAL(X)
Definition: Random.h:29
long Random::integer ( void  )
inline

Definition at line 214 of file Random.h.

References skip(), and UINT64_LITERAL.

Referenced by reorder().

214  {
215  skip();
216  return ( ( rand48_seed >> 17 ) & UINT64_LITERAL(0x000000007fffffff) );
217  }
void skip(void)
Definition: Random.h:72
#define UINT64_LITERAL(X)
Definition: Random.h:29
template<class Elem >
void Random::reorder ( Elem *  a,
int  n 
)
inline

Definition at line 220 of file Random.h.

References integer().

Referenced by WorkDistrib::buildNodeAwarePeOrdering(), ComputePmeMgr::initialize(), ComputePmeMgr::initialize_pencils(), PmePencil< CBase_PmeZPencil >::order_init(), HomePatch::registerProxy(), and ParallelIOMgr::sendAtomsToHomePatchProcs().

220  {
221  for ( int i = 0; i < (n-1); ++i ) {
222  int ie = i + ( integer() % (n-i) );
223  if ( ie == i ) continue;
224  const Elem e = a[ie];
225  a[ie] = a[i];
226  a[i] = e;
227  }
228  }
long integer(void)
Definition: Random.h:214
void Random::skip ( void  )
inline

Definition at line 72 of file Random.h.

References RAND48_MASK.

Referenced by integer(), split(), and uniform().

72  {
73  rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) & RAND48_MASK;
74  }
#define RAND48_MASK
Definition: Random.h:35
void Random::split ( int  iStream,
int  numStreams 
)
inline

Definition at line 77 of file Random.h.

References skip().

Referenced by Controller::Controller(), Sequencer::Sequencer(), and Node::startup().

77  {
78 
79  int i;
80 
81  // make sure that numStreams is odd to ensure maximum period
82  numStreams |= 1;
83 
84  // iterate to get to the correct stream
85  for ( i = 0; i < iStream; ++i ) skip();
86 
87  // save seed and add so we can use skip() for our calculations
88  uint64_t save_seed = rand48_seed;
89 
90  // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
91  rand48_seed = rand48_add;
92  for ( i = 1; i < numStreams; ++i ) skip();
93  uint64_t new_add = rand48_seed;
94 
95  // calculate a = a^numStreams
96  rand48_seed = rand48_mult;
97  rand48_add = 0;
98  for ( i = 1; i < numStreams; ++i ) skip();
99  rand48_mult = rand48_seed;
100 
101  rand48_add = new_add;
102  rand48_seed = save_seed;
103 
104  second_gaussian = 0;
105  second_gaussian_waiting = 0;
106  }
void skip(void)
Definition: Random.h:72
BigReal Random::sum_of_squared_gaussians ( int64_t  num_gaussians)
inline

Return the sum of i.i.d. squared Gaussians.

The sum of k squared standard normal random variables is equivalent to drawing a value from a chi-squared distribution with the shape given as the number of Gaussians. That is,

X ~ sum_n=1^k N(0, 1)^2 ~ chi^2(k).

This is in turn a special case of the Gamma distribution with shape k/2 and scale 2 (we could also use rate = 1/scale)

X ~ chi^2(k) ~ Gamma(k/2, 2) ~ 2*Gamma(k/2, 1).

The second relation follows from the scaling property of the Gamma distribution. Furthermore, when a Gamma distribution has unit scale and a large shape, it can be well-approximated by a normal distribution with mean and variance equal to the shape. Thus,

X ~ 2*Gamma(k/2, 1) ~= 2*N(k/2, sqrt(k/2)).

A quick numerical test shows that the mean integrated square error for this approximation is <10^-5 for shape >30 and <10^-6 for shape >100. We'll be conservative and use the latter cutoff.

We thus have three cases for k Gaussians:

0 < k <= 2 - just brute force generate and sum the Gaussians 2 < k <= 200 - use a (slightly modified) version of the Gamma distribution algorithm from Marsaglia and Tsang that is implemented in the GNU Science Library (GSL) else - use a single Gaussian distribution

The brute force method is almost certainly the slowest method, even for k = 3. The rigorous method takes about 150% as long as the approximate method (but this is invariant to the value of k).

Parameters
num_gaussiansA positive integer number of Gaussians
Returns
a random variable equal to the sum of num_gaussians squared standard normal variables

Definition at line 178 of file Random.h.

References gaussian(), uniform(), and z.

Referenced by Controller::stochRescaleCoefficient().

178  {
179  BigReal z, u, v;
180 
181  if (num_gaussians <= 2) {
182  v = 0.0;
183  for(int i=0; i<num_gaussians; ++i) {
184  z = gaussian();
185  v += z*z;
186  }
187  return v;
188  } else if (2 < num_gaussians && num_gaussians <= 200) {
189  const BigReal d = 0.5*num_gaussians - 1./3;
190  const BigReal c = 1. / (3*sqrt(d));
191  const BigReal zmin = -1. / c;
192  do {
193  do {
194  z = gaussian();
195  }
196  while(z <= zmin);
197  u = uniform();
198  v = (1 + c*z); v = v*v*v; // v = (1 + c*z)^3
199  }
200  while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
201  return 2*d*v;
202  } else {
203  return num_gaussians + sqrt(2*num_gaussians)*gaussian();
204  }
205  }
BigReal uniform(void)
Definition: Random.h:109
BigReal gaussian(void)
Definition: Random.h:116
gridSize z
double BigReal
Definition: common.h:112
BigReal Random::uniform ( void  )
inline

Definition at line 109 of file Random.h.

References skip(), and UINT64_LITERAL.

Referenced by Controller::adaptTempUpdate(), gaussian(), sum_of_squared_gaussians(), and uniform_array_f().

109  {
110  skip();
111  const double exp48 = ( 1.0 / (double)(UINT64_LITERAL(1) << 48) );
112  return ( (double) rand48_seed * exp48 );
113  }
void skip(void)
Definition: Random.h:72
#define UINT64_LITERAL(X)
Definition: Random.h:29
void Random::uniform_array_f ( float *  a,
int  n 
)
inline

Fill length n array from uniform distribution.

Definition at line 233 of file Random.h.

References uniform().

233  {
234  // for now just call uniform()
235  // ultimately we will provide faster implementation
236  for (int i=0; i < n; i++) {
237  a[i] = (float)uniform();
238  }
239  }
BigReal uniform(void)
Definition: Random.h:109

The documentation for this class was generated from the following file: