Random Class Reference

#include <Random.h>

List of all members.

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 (int num_gaussians)
 Return the sum of i.i.d. squared Gaussians.
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.

00050                {
00051     init(0);
00052     rand48_seed = RAND48_SEED;
00053   }

Random::Random ( unsigned long  seed  )  [inline]

Definition at line 56 of file Random.h.

References init().

00056                              {
00057     init(seed);
00058   }


Member Function Documentation

BigReal Random::gaussian ( void   )  [inline]

Definition at line 116 of file Random.h.

References uniform().

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

00116                          {
00117     BigReal fac, r, v1, v2;
00118 
00119     if (second_gaussian_waiting) {
00120       second_gaussian_waiting = 0;
00121       return second_gaussian;
00122     } else {
00123       r = 2.;                 // r >= 1.523e-8 ensures abs result < 6
00124       while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
00125         v1 = 2.0 * uniform() - 1.0;
00126         v2 = 2.0 * uniform() - 1.0;
00127         r = v1*v1 + v2*v2;
00128       }
00129       fac = sqrt(-2.0 * log(r)/r);
00130       // now make the Box-Muller transformation to get two normally
00131       // distributed random numbers. Save one and return the other.
00132       second_gaussian_waiting = 1;
00133       second_gaussian = v1 * fac;
00134       return v2 * fac;
00135     }
00136   }

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().

00244                                          {
00245     // for now just call gaussian()
00246     // ultimately we will provide faster implementation
00247     for (int i=0;  i < n;  i++) {
00248       a[i] = (float)gaussian();
00249     }
00250   }

Vector Random::gaussian_vector ( void   )  [inline]

Definition at line 208 of file Random.h.

References gaussian().

Referenced by Controller::langevinPiston1(), Controller::langevinPiston2(), Sequencer::langevinVelocities(), Sequencer::langevinVelocitiesBBK2(), Sequencer::reassignVelocities(), and Sequencer::reinitVelocities().

00208                                {
00209     return Vector( gaussian(), gaussian(), gaussian() );
00210   }

void Random::init ( unsigned long  seed  )  [inline]

Definition at line 61 of file Random.h.

References INT64_LITERAL, RAND48_ADD, RAND48_MULT, and RAND48_SEED.

Referenced by Random().

00061                                 {
00062     second_gaussian = 0;
00063     second_gaussian_waiting = 0;
00064     rand48_seed = seed & INT64_LITERAL(0x00000000ffffffff);
00065     rand48_seed = rand48_seed << 16;
00066     rand48_seed |= RAND48_SEED & INT64_LITERAL(0x0000ffff);
00067     rand48_mult = RAND48_MULT;
00068     rand48_add = RAND48_ADD;
00069   }

long Random::integer ( void   )  [inline]

Definition at line 214 of file Random.h.

References INT64_LITERAL, and skip().

Referenced by reorder().

00214                      {
00215     skip();
00216     return ( ( rand48_seed >> 17 ) & INT64_LITERAL(0x000000007fffffff) );
00217   }

template<class Elem>
void Random::reorder ( Elem *  a,
int  n 
) [inline]

Definition at line 220 of file Random.h.

References integer().

Referenced by ComputePmeMgr::initialize(), and OptPmePencil< CBase_OptPmeXPencil >::order_init().

00220                                                      {
00221     for ( int i = 0; i < (n-1); ++i ) {
00222       int ie = i + ( integer() % (n-i) );
00223       if ( ie == i ) continue;
00224       const Elem e = a[ie];
00225       a[ie] = a[i];
00226       a[i] = e;
00227     }
00228   }

void Random::skip ( void   )  [inline]

Definition at line 72 of file Random.h.

References RAND48_MASK.

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

00072                   {
00073     rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) & RAND48_MASK;
00074   }

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().

00077                                           {
00078 
00079     int i;
00080 
00081     // make sure that numStreams is odd to ensure maximum period
00082     numStreams |= 1;
00083 
00084     // iterate to get to the correct stream
00085     for ( i = 0; i < iStream; ++i ) skip();
00086 
00087     // save seed and add so we can use skip() for our calculations
00088     int64 save_seed = rand48_seed;
00089 
00090     // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
00091     rand48_seed = rand48_add;
00092     for ( i = 1; i < numStreams; ++i ) skip();
00093     int64 new_add = rand48_seed;
00094 
00095     // calculate a = a^numStreams
00096     rand48_seed = rand48_mult;
00097     rand48_add  = 0;
00098     for ( i = 1; i < numStreams; ++i ) skip();
00099     rand48_mult = rand48_seed;
00100 
00101     rand48_add  = new_add;
00102     rand48_seed = save_seed;
00103 
00104     second_gaussian = 0;
00105     second_gaussian_waiting = 0;
00106   }

BigReal Random::sum_of_squared_gaussians ( int  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_gaussians A 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::stochRescaleVelocities().

00178                                                       {
00179     BigReal z, u, v;
00180 
00181     if (num_gaussians <= 2) {
00182         v = 0.0;
00183         for(int i=0; i<num_gaussians; ++i) {
00184           z = gaussian();
00185           v += z*z;
00186         }
00187         return v;
00188     } else if (2 < num_gaussians && num_gaussians <= 200) {
00189       const BigReal d = 0.5*num_gaussians - 1./3;
00190       const BigReal c = 1. / (3*sqrt(d));
00191       const BigReal zmin = -1. / c;
00192       do {
00193         do {
00194           z = gaussian();
00195         }
00196         while(z <= zmin);
00197         u = uniform();
00198         v = (1 + c*z); v = v*v*v; // v = (1 + c*z)^3
00199       }
00200       while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
00201       return 2*d*v;
00202     } else {
00203       return num_gaussians + sqrt(2*num_gaussians)*gaussian();
00204     }
00205   }

BigReal Random::uniform ( void   )  [inline]

Definition at line 109 of file Random.h.

References INT64_LITERAL, and skip().

Referenced by gaussian(), SELF(), sum_of_squared_gaussians(), and uniform_array_f().

00109                         {
00110     skip();
00111     const double exp48 = ( 1.0 / (double)(INT64_LITERAL(1) << 48) );
00112     return ( (double) rand48_seed * exp48 );
00113   }

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().

00233                                         {
00234     // for now just call uniform()
00235     // ultimately we will provide faster implementation
00236     for (int i=0;  i < n;  i++) {
00237       a[i] = (float)uniform();
00238     }
00239   }


The documentation for this class was generated from the following file:
Generated on Sat Jun 23 01:17:22 2018 for NAMD by  doxygen 1.4.7