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() [1/2]

Random::Random ( void  )

inline

Definition at line 50 of file Random.h.

References init(), and RAND48_SEED.

               {
    init(0);
    rand48_seed = RAND48_SEED;
  }

Random() [2/2]

Random::Random ( unsigned long  seed )

inline

Definition at line 56 of file Random.h.

References init().

                             {
    init(seed);
  }

Member Function Documentation

gaussian()

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

                         {
#ifdef DEBUG_RANDOM
    static int count = 0;
#endif
    BigReal fac, r, v1, v2;

    if (second_gaussian_waiting) {
      second_gaussian_waiting = 0;
#ifdef DEBUG_RANDOM
      count++;
      fprintf(stderr, "RANDOM  %6d  %12.9f\n", count, second_gaussian);
#endif
      return second_gaussian;
    } else {
      r = 2.;                 // r >= 1.523e-8 ensures abs result < 6
      while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
        v1 = 2.0 * uniform() - 1.0;
        v2 = 2.0 * uniform() - 1.0;
        r = v1*v1 + v2*v2;
      }
      fac = sqrt(-2.0 * log(r)/r);
      // now make the Box-Muller transformation to get two normally
      // distributed random numbers. Save one and return the other.
      second_gaussian_waiting = 1;
      second_gaussian = v1 * fac;
#ifdef DEBUG_RANDOM
      count++;
      fprintf(stderr, "RANDOM  %6d  %12.9f\n", count, v2 * fac);
#endif
      return v2 * fac;
    }
  }

gaussian_array_f()

void Random::gaussian_array_f ( float *  a, int  n  )

inline

Fill length n array from standard Gaussian distribution.

Definition at line 255 of file Random.h.

References gaussian().

Referenced by Sequencer::langevinVelocitiesBBK2_SOA().

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

gaussian_vector()

Vector Random::gaussian_vector ( void  )

inline

Definition at line 219 of file Random.h.

References gaussian().

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

                               {
    return Vector( gaussian(), gaussian(), gaussian() );
  }

init()

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

                                {
    second_gaussian = 0;
    second_gaussian_waiting = 0;
    rand48_seed = seed & UINT64_LITERAL(0x00000000ffffffff);
    rand48_seed = rand48_seed << 16;
    rand48_seed |= RAND48_SEED & UINT64_LITERAL(0x0000ffff);
    rand48_mult = RAND48_MULT;
    rand48_add = RAND48_ADD;
  }

integer()

long Random::integer ( void  )

inline

Definition at line 225 of file Random.h.

References skip(), and UINT64_LITERAL.

Referenced by reorder().

                     {
    skip();
    return ( ( rand48_seed >> 17 ) & UINT64_LITERAL(0x000000007fffffff) );
  }

reorder()

template<class Elem >

void Random::reorder ( Elem *  a, int  n  )

inline

Definition at line 231 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().

                                                     {
    for ( int i = 0; i < (n-1); ++i ) {
      int ie = i + ( integer() % (n-i) );
      if ( ie == i ) continue;
      const Elem e = a[ie];
      a[ie] = a[i];
      a[i] = e;
    }
  }

skip()

void Random::skip ( void  )

inline

Definition at line 72 of file Random.h.

References RAND48_MASK.

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

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

split()

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

                                          {

    int i;

    // make sure that numStreams is odd to ensure maximum period
    numStreams |= 1;

    // iterate to get to the correct stream
    for ( i = 0; i < iStream; ++i ) skip();

    // save seed and add so we can use skip() for our calculations
    uint64_t save_seed = rand48_seed;

    // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
    rand48_seed = rand48_add;
    for ( i = 1; i < numStreams; ++i ) skip();
    uint64_t new_add = rand48_seed;

    // calculate a = a^numStreams
    rand48_seed = rand48_mult;
    rand48_add  = 0;
    for ( i = 1; i < numStreams; ++i ) skip();
    rand48_mult = rand48_seed;

    rand48_add  = new_add;
    rand48_seed = save_seed;

    second_gaussian = 0;
    second_gaussian_waiting = 0;
  }

sum_of_squared_gaussians()

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_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 189 of file Random.h.

References gaussian(), and uniform().

Referenced by Controller::stochRescaleCoefficient().

                                                          {
    BigReal z, u, v;

    if (num_gaussians <= 2) {
        v = 0.0;
        for(int i=0; i<num_gaussians; ++i) {
          z = gaussian();
          v += z*z;
        }
        return v;
    } else if (2 < num_gaussians && num_gaussians <= 200) {
      const BigReal d = 0.5*num_gaussians - 1./3;
      const BigReal c = 1. / (3*sqrt(d));
      const BigReal zmin = -1. / c;
      do {
        do {
          z = gaussian();
        }
        while(z <= zmin);
        u = uniform();
        v = (1 + c*z); v = v*v*v; // v = (1 + c*z)^3
      }
      while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
      return 2*d*v;
    } else {
      return num_gaussians + sqrt(2*num_gaussians)*gaussian();
    }
  }

uniform()

BigReal Random::uniform ( void  )

inline

Definition at line 109 of file Random.h.

References skip(), and UINT64_LITERAL.

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

                        {
    skip();
    const double exp48 = ( 1.0 / (double)(UINT64_LITERAL(1) << 48) );
    return ( (double) rand48_seed * exp48 );
  }

uniform_array_f()

void Random::uniform_array_f ( float *  a, int  n  )

inline

Fill length n array from uniform distribution.

Definition at line 244 of file Random.h.

References uniform().

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

---

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

• Random.h