/home/fwph/code/wurde/rde/core/Random.H

Go to the documentation of this file.

#ifndef Random_H
#define Random_H

#include <cmath>
#include <unistd.h>
#include <fcntl.h>




class MTRand {
public:
  MTRand(const unsigned long &oneSeed) {setSeed(oneSeed);}


  void setSeed(unsigned long seed);

  void randomize();
  


  virtual double mean() const = 0;
  virtual double variance() const = 0;
  virtual double stdDev() const = 0;

private:
  typedef unsigned long uint32;
  
protected:
  double rand() {return double(randInt()) * 2.3283064370807974e-10;}

  uint32 randInt();
  uint32 randInt(const uint32 &n);

private:
  enum {N = 624, M = 397, MAGIC = 0x9908b0dfU};
  
  uint32 state[N];
  uint32 *pNext;
  int left;
  
  void reload();

  uint32 hiBit(const uint32 &u) const {return u & 0x80000000U;}
  uint32 loBit(const uint32 &u) const {return u & 0x00000001U;}
  uint32 loBits(const uint32 &u) const {return u & 0x7fffffffU;}
  uint32 mixBits(const uint32 &u, const uint32 &v) const {
    return hiBit(u) | loBits(v);}
  uint32 twist(const uint32 &m, const uint32 &s0, const uint32 &s1) const {
    return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? MAGIC : 0U);}

  // Not allowed to use these
  MTRand(const MTRand &mtrand);
  MTRand &operator=(const MTRand &mtrand);
};



class RealRandom : public MTRand {
public:
  RealRandom(const unsigned long seed) :MTRand(seed) {}




  RealRandom &operator>>(double &random_number) {
    random_number = number(); return *this;}

  double operator()() {return number();}

  operator double() {return number();}


  virtual double number() = 0;


  virtual double mean() const = 0;
  virtual double variance() const = 0;
  virtual double stdDev() const = 0;
};



class IntegerRandom : public MTRand {
public:
  IntegerRandom(const unsigned long seed) :MTRand(seed) {}




  IntegerRandom &operator>>(int &random_number) {
    random_number = number(); return *this;}

  int operator()() {return number();}

  operator int() {return number();}


  virtual int number() = 0;


  virtual double mean() const = 0;
  virtual double variance() const = 0;
  virtual double stdDev() const = 0;
};



class UniformRandomSimple : public RealRandom {
public:
  UniformRandomSimple(const unsigned long seed = 1UL) :RealRandom(seed) {}


  double number() {return MTRand::rand();}


  double mean() const {return 0.5;}
  double variance() const {return (1.0 / 12.0);}
  double stdDev() const {return sqrt(1.0 / 12.0);}
};



class UniformRandom : public RealRandom {
public:
  UniformRandom(const double low = 1.0, const double high = 0.0,
                const unsigned long seed = 1UL)
    :RealRandom(seed), m_low((low < high)?low:high),
     m_range((low < high)?(high - low):(low - high)) {}


  double number() {return (MTRand::rand() * m_range) + m_low;}


  double mean() const {return ((m_range / 2.0) + m_low);}
  double variance() const {return ((m_range * m_range) / 12.0);}
  double stdDev() const {return sqrt(variance());}



  double lowValue() const {return m_low;}
  double highValue() const {return (m_low + m_range);}
  double range() const {return m_range;}



  void setLowValue(const double low) {m_low = low;}
  void setHighValue(const double high) {m_range = high - m_low;}  
  void setValues(const double low, const double high) {
    setLowValue(low); setHighValue(high);}

private:
  double m_low;
  double m_range;
};



class NormalRandom : public RealRandom {
public:
  NormalRandom(const double mean = 0.0, const double variance = 1.0,
               const unsigned long seed = 1UL)
    :RealRandom(seed), m_mean(mean), m_std_dev(sqrt(fabs(variance))) {}
  

  double number();


  double mean() const {return m_mean;}
  double variance() const {return (m_std_dev * m_std_dev);}
  double stdDev() const {return m_std_dev;}



  void setMean(const double mean) {m_mean = mean;}
  void setVariance(const double variance) {m_std_dev = sqrt(fabs(variance));}
  void setStdDev(const double std_dev) {m_std_dev = fabs(std_dev);}
  void setParameters(const double mean, const double variance) {
      setMean(mean); setVariance(variance);}

private:
  double m_mean;
  double m_std_dev;
};



class ExponentialRandom : public RealRandom {
public:
  ExponentialRandom(const double lambda = 1.0, const unsigned long seed = 1UL)
    :RealRandom(seed), m_lambda(lambda) {}
  

  double number();


  double mean() const {return (1.0 / m_lambda);}
  double variance() const {return (1.0 / (m_lambda * m_lambda));}
  double stdDev() const {return (1.0 / m_lambda);}



  double lambda() const {return m_lambda;}



  void setMean(const double mean) {m_lambda = (1.0 / mean);}
  void setVariance(const double variance) {m_lambda = sqrt(1.0 / variance);}
  void setStdDev(const double std_dev) {m_lambda = (1.0 / std_dev);}
  void setLambda(const double lambda) {m_lambda = lambda;}
  void setParameters(const double mean, const double variance) {
      setMean(mean); setVariance(variance);}

private:
  double m_lambda;
};



class DiscreteUniformRandom : public IntegerRandom {
public:
  DiscreteUniformRandom(const int low, const int high,
                       const unsigned long seed = 1UL)
    :IntegerRandom(seed), m_low((low < high)?low:high),
     m_range((unsigned long)((low < high)?(high - low):(low - high))) {}
  DiscreteUniformRandom(const int range, const unsigned long seed = 1UL)
    :IntegerRandom(seed), m_low(0), m_range(range - 1) {}


  int number() {return MTRand::randInt(m_range) + m_low;}


  double mean() const {return (((double)m_range / 2.0));}
  double variance() const {return ((double)(m_range * m_range) / 12.0);}
  double stdDev() const {return sqrt(variance());}



  int lowValue() const {return m_low;}
  int highValue() const {return (m_low + m_range);}
  int range() const {return m_range + 1;}



  void setLowValue(const int low) {m_low = low;}
  void setHighValue(const int high) {m_range = (unsigned long)(high - m_low);}
  void setValues(const int low, const int high) {
    setLowValue(low); setHighValue(high);}

private:
  int m_low;
  unsigned long m_range;
};



class BernoulliRandom : public IntegerRandom {
public:
  BernoulliRandom(const double success = 0.5, const unsigned long seed = 1UL)
    :IntegerRandom(seed), m_success(success) {}


  int number() {return (rand() <= m_success)?1:0;}


  double mean() const {return m_success;}
  double variance() const {return (m_success * (1.0 - m_success));}
  double stdDev() const {return sqrt(variance());}



  double success() const {return m_success;}
  double probSuccess() const {return success();}
  double failure() const {return (1.0 - m_success);}
  double probFailure() const {return failure();}



  void setProbSuccess(const double success) {m_success = success;}
  void setProbFailure(const double failure) {m_success = 1.0 - failure;}

private:
  double m_success;
};


// This is ugly, but half of the cost of these calls is in the function
// call overhead.


inline MTRand::uint32 MTRand::randInt() {
  if(left == 0)
    reload();
  --left;
                
  register uint32 s1;
  s1 = *pNext++;
  s1 ^= (s1 >> 11);
  s1 ^= (s1 <<  7) & 0x9d2c5680U;
  s1 ^= (s1 << 15) & 0xefc60000U;
  return (s1 ^ (s1 >> 18));
}


inline MTRand::uint32 MTRand::randInt(const uint32 &n) {
        // RLG: this is a problem; ~0 has the sign bit set.
  // Find which bits are used in n
        //  uint32 used = ~0;
        long int used = ~0;
  for(uint32 m = n; m; used <<= 1, m >>= 1 ) {}
  used = ~used;
        
  // Draw numbers until one is found in [0,n]
  uint32 i;
  do
    i = randInt() & used;  // toss unused bits to shorten search
  while(i > n);
  return i;
}


inline void MTRand::setSeed(unsigned long seed) {
  // Seed the generator with a simple uint32
  register uint32 *s;
  register int i;
  for(i = N, s = state; i--;
      *s = seed & 0xffff0000,
        *s++ |= ((seed *= 69069U)++ & 0xffff0000) >> 16,
        (seed *= 69069U)++) {}  // hard to read, but fast
  reload();
}


inline void MTRand::randomize() {
  unsigned long seed_number;

  int fp = open("/dev/urandom", O_RDONLY);
  read(fp, &seed_number, sizeof(seed_number));
  close(fp);

  MTRand::setSeed(seed_number);
}



inline void MTRand::reload() {
  // Generate N new values in state
  // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
  register uint32 *p = state;
  register int i;
  for(i = N - M; i--;) {
                //              *p++ = twist(p[M], p[0], p[1]);    // RLG: Precedence poorly defined
                *p = twist(p[M], p[0], p[1]);
                p++;
        }
  for(i = M; --i;) {
                //              *p++ = twist(p[M-N], p[0], p[1]);  // RLG: Precedence poorly defined.
                *p = twist(p[M-N], p[0], p[1]);
                p++;
        }
  *p = twist(p[M-N], p[0], state[0]);
  left = N, pNext = state;
}


// We could cache another number m_std_dev * r * sin(theta) + m_mean, but
// this introduces non-randomness into the stream of numbers, although it
// speeds things up a bit.
inline double NormalRandom::number() {
  double theta = 2.0 * M_PI * MTRand::rand();
  double r = sqrt(-2.0 * log(MTRand::rand()));
  
  return ((m_std_dev * r * cos(theta)) + m_mean);
}


#endif

---