--- a/src/core/random-variable.cc Wed Feb 04 17:47:29 2009 -0800 +++ a/src/core/random-variable.cc Wed May 27 15:48:57 2009 +0200 @@ -16,6 +16,7 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // Author: Rajib Bhattacharjea +// Author: Hadi Arbabi // #include @@ -32,6 +33,8 @@ #include "assert.h" +#include "config.h" +#include "integer.h" #include "random-variable.h" #include "rng-stream.h" #include "fatal-error.h" @@ -41,6 +44,44 @@ namespace ns3{ //----------------------------------------------------------------------------- +// Seed Manager +//----------------------------------------------------------------------------- + +uint32_t SeedManager::GetSeed() +{ + uint32_t s[6]; + RngStream::GetPackageSeed (s); + NS_ASSERT( + s[0] == s[1] && + s[0] == s[2] && + s[0] == s[3] && + s[0] == s[4] && + s[0] == s[5] + ); + return s[0]; +} + +void SeedManager::SetSeed(uint32_t seed) +{ + Config::SetGlobal("RngSeed", IntegerValue(seed)); +} + +void SeedManager::SetRun(uint32_t run) +{ + Config::SetGlobal("RngRun", IntegerValue(run)); +} + +uint32_t SeedManager::GetRun() +{ + return RngStream::GetPackageRun (); +} + +bool SeedManager::CheckSeed (uint32_t seed) +{ + return RngStream::CheckSeed(seed); +} + +//----------------------------------------------------------------------------- //----------------------------------------------------------------------------- // RandomVariableBase methods @@ -54,72 +95,23 @@ virtual double GetValue() = 0; virtual uint32_t GetInteger(); virtual RandomVariableBase* Copy(void) const = 0; - virtual void GetSeed(uint32_t seed[6]); - static void UseDevRandom(bool udr = true); - static void UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, - uint32_t s3, uint32_t s4, uint32_t s5); - static void SetRunNumber(uint32_t n); -private: - static void GetRandomSeeds(uint32_t seeds[6]); -private: - static bool useDevRandom; // True if using /dev/random desired - static bool globalSeedSet; // True if global seed has been specified - static int devRandom; // File handle for /dev/random - static uint32_t globalSeed[6]; // The global seed to use - friend class RandomVariableInitializer; protected: - static unsigned long heuristic_sequence; - static RngStream* m_static_generator; - static uint32_t runNumber; - static void Initialize(); // Initialize the RNG system - static bool initialized; // True if package seed is set RngStream* m_generator; //underlying generator being wrapped }; - - -bool RandomVariableBase::initialized = false; // True if RngStream seed set -bool RandomVariableBase::useDevRandom = false; // True if use /dev/random -bool RandomVariableBase::globalSeedSet = false; // True if GlobalSeed called -int RandomVariableBase::devRandom = -1; -uint32_t RandomVariableBase::globalSeed[6]; -unsigned long RandomVariableBase::heuristic_sequence; -RngStream* RandomVariableBase::m_static_generator = 0; -uint32_t RandomVariableBase::runNumber = 0; - -//the static object random_variable_initializer initializes the static members -//of RandomVariable -static class RandomVariableInitializer -{ - public: - RandomVariableInitializer() - { -// RandomVariableBase::Initialize(); // sets the static package seed -// RandomVariableBase::m_static_generator = new RngStream(); -// RandomVariableBase::m_static_generator->InitializeStream(); - } - ~RandomVariableInitializer() - { - delete RandomVariableBase::m_static_generator; - } -} random_variable_initializer; - RandomVariableBase::RandomVariableBase() : m_generator(NULL) { -// m_generator = new RngStream(); -// m_generator->InitializeStream(); -// m_generator->ResetNthSubstream(RandomVariableBase::runNumber); } RandomVariableBase::RandomVariableBase(const RandomVariableBase& r) :m_generator(0) { - if(r.m_generator) - { - m_generator = new RngStream(*r.m_generator); - } + if (r.m_generator) + { + m_generator = new RngStream(*r.m_generator); + } } RandomVariableBase::~RandomVariableBase() @@ -132,115 +124,7 @@ return (uint32_t)GetValue(); } -void RandomVariableBase::UseDevRandom(bool udr) -{ - RandomVariableBase::useDevRandom = udr; -} - -void RandomVariableBase::GetSeed(uint32_t seed[6]) -{ - if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } - m_generator->GetState(seed); -} - -//----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- -// RandomVariableBase static methods -void RandomVariableBase::UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, - uint32_t s3, uint32_t s4, uint32_t s5) -{ - if (RandomVariableBase::globalSeedSet) - { - NS_FATAL_ERROR ("Random number generator already initialized! " - "Call to RandomVariableBase::UseGlobalSeed() ignored"); - } - RandomVariableBase::globalSeed[0] = s0; - RandomVariableBase::globalSeed[1] = s1; - RandomVariableBase::globalSeed[2] = s2; - RandomVariableBase::globalSeed[3] = s3; - RandomVariableBase::globalSeed[4] = s4; - RandomVariableBase::globalSeed[5] = s5; - if (!RngStream::CheckSeed(RandomVariableBase::globalSeed)) - NS_FATAL_ERROR("Invalid seed"); - - RandomVariableBase::globalSeedSet = true; -} - -void RandomVariableBase::Initialize() -{ - if (RandomVariableBase::initialized) return; // Already initialized and seeded - RandomVariableBase::initialized = true; - if (!RandomVariableBase::globalSeedSet) - { // No global seed, try a random one - GetRandomSeeds(globalSeed); - } - // Seed the RngStream package - RngStream::SetPackageSeed(globalSeed); -} - -void RandomVariableBase::GetRandomSeeds(uint32_t seeds[6]) -{ - // Check if /dev/random exists - if (RandomVariableBase::useDevRandom && RandomVariableBase::devRandom < 0) - { - RandomVariableBase::devRandom = open("/dev/random", O_RDONLY); - } - if (RandomVariableBase::devRandom > 0) - { // Use /dev/random - while(true) - { - for (int i = 0; i < 6; ++i) - { - ssize_t bytes_read = read (RandomVariableBase::devRandom, - &seeds[i], sizeof (seeds[i])); - if (bytes_read != sizeof (seeds[i])) - { - NS_FATAL_ERROR ("Read from /dev/random failed"); - } - } - if (RngStream::CheckSeed(seeds)) break; // Got a valid one - } - } - else - { // Seed from time of day (code borrowed from ns2 random seeding) - // Thanks to John Heidemann for this technique - while(true) - { - timeval tv; - gettimeofday(&tv, 0); - seeds[0] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - gettimeofday(&tv, 0); - seeds[1] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - gettimeofday(&tv, 0); - seeds[2] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - gettimeofday(&tv, 0); - seeds[3] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - gettimeofday(&tv, 0); - seeds[4] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - gettimeofday(&tv, 0); - seeds[5] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8)) - & 0x7fffffff; - if (RngStream::CheckSeed(seeds)) break; // Got a valid one - } - } -} - -void RandomVariableBase::SetRunNumber(uint32_t n) -{ - runNumber = n; -} - - +//------------------------------------------------------- RandomVariable::RandomVariable() : m_variable (0) @@ -277,33 +161,14 @@ { return m_variable->GetInteger (); } -void -RandomVariable::GetSeed(uint32_t seed[6]) const -{ - return m_variable->GetSeed (seed); -} -void -RandomVariable::UseDevRandom(bool udr) -{ - RandomVariableBase::UseDevRandom (udr); -} -void -RandomVariable::UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, - uint32_t s3, uint32_t s4, uint32_t s5) -{ - RandomVariableBase::UseGlobalSeed (s0, s1, s2, s3, s4, s5); -} -void -RandomVariable::SetRunNumber(uint32_t n) -{ - RandomVariableBase::SetRunNumber (n); -} + RandomVariableBase * RandomVariable::Peek (void) const { return m_variable; } + ATTRIBUTE_VALUE_IMPLEMENT (RandomVariable); ATTRIBUTE_CHECKER_IMPLEMENT (RandomVariable); @@ -335,15 +200,14 @@ * \return A value between low and high values specified by the constructor */ virtual double GetValue(); + + /** + * \return A value between low and high values specified by parameters + */ + virtual double GetValue(double s, double l); + virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param s Low end of the range - * \param l High end of the range - * \return A uniformly distributed random number between s and l - */ - static double GetSingleValue(double s, double l); private: double m_min; double m_max; @@ -372,47 +236,49 @@ double UniformVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } return m_min + m_generator->RandU01() * (m_max - m_min); } +double UniformVariableImpl::GetValue(double s, double l) +{ + if(!m_generator) + { + m_generator = new RngStream(); + } + return s + m_generator->RandU01() * (l-s); +} + RandomVariableBase* UniformVariableImpl::Copy() const { return new UniformVariableImpl(*this); } -double UniformVariableImpl::GetSingleValue(double s, double l) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - return s + m_static_generator->RandU01() * (l - s);; -} - UniformVariable::UniformVariable() : RandomVariable (UniformVariableImpl ()) {} UniformVariable::UniformVariable(double s, double l) : RandomVariable (UniformVariableImpl (s, l)) {} -double -UniformVariable::GetSingleValue(double s, double l) + +double UniformVariable::GetValue(void) const { - return UniformVariableImpl::GetSingleValue (s, l); + return this->RandomVariable::GetValue (); } +double UniformVariable::GetValue(double s, double l) +{ + return ((UniformVariableImpl*)Peek())->GetValue(s,l); +} + +uint32_t UniformVariable::GetInteger (uint32_t s, uint32_t l) +{ + NS_ASSERT(s <= l); + return static_cast( GetValue(s, l+1) ); +} //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -623,13 +489,7 @@ */ virtual double GetValue(); virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param m The mean of the distribution from which the return value is drawn - * \param b The upper bound value desired, beyond which values get clipped - * \return A random number from an exponential distribution with mean m - */ - static double GetSingleValue(double m, double b=0); + private: double m_mean; // Mean value of RV double m_bound; // Upper bound on value (if non-zero) @@ -649,37 +509,22 @@ double ExponentialVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } - double r = -m_mean*log(m_generator->RandU01()); - if (m_bound != 0 && r > m_bound) return m_bound; - return r; + { + m_generator = new RngStream(); + } + while(1) + { + double r = -m_mean*log(m_generator->RandU01()); + if (m_bound == 0 || r <= m_bound) return r; + //otherwise, try again + } } RandomVariableBase* ExponentialVariableImpl::Copy() const { return new ExponentialVariableImpl(*this); } -double ExponentialVariableImpl::GetSingleValue(double m, double b/*=0*/) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - double r = -m*log(m_static_generator->RandU01()); - if (b != 0 && r > b) return b; - return r; -} ExponentialVariable::ExponentialVariable() : RandomVariable (ExponentialVariableImpl ()) @@ -690,11 +535,6 @@ ExponentialVariable::ExponentialVariable(double m, double b) : RandomVariable (ExponentialVariableImpl (m, b)) {} -double -ExponentialVariable::GetSingleValue(double m, double b) -{ - return ExponentialVariableImpl::GetSingleValue (m, b); -} //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -743,17 +583,7 @@ */ virtual double GetValue(); virtual RandomVariableBase* Copy() const; -public: - /** - * \param m The mean value of the distribution from which the return value - * is drawn. - * \param s The shape parameter of the distribution from which the return - * value is drawn. - * \param b The upper bound to which to restrict return values - * \return A random number from a Pareto distribution with mean m and shape - * parameter s. - */ - static double GetSingleValue(double m, double s, double b=0); + private: double m_mean; // Mean value of RV double m_shape; // Shape parameter @@ -778,20 +608,17 @@ double ParetoVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } double scale = m_mean * ( m_shape - 1.0) / m_shape; - double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape))); - if (m_bound != 0 && r > m_bound) return m_bound; - return r; + while(1) + { + double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape))); + if (m_bound == 0 || r <= m_bound) return r; + //otherwise, try again + } } RandomVariableBase* ParetoVariableImpl::Copy() const @@ -799,20 +626,6 @@ return new ParetoVariableImpl(*this); } -double ParetoVariableImpl::GetSingleValue(double m, double s, double b/*=0*/) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - double scale = m * ( s - 1.0) / s; - double r = (scale * ( 1.0 / pow(m_static_generator->RandU01(), 1.0 / s))); - if (b != 0 && r > b) return b; - return r; -} - ParetoVariable::ParetoVariable () : RandomVariable (ParetoVariableImpl ()) {} @@ -825,11 +638,6 @@ ParetoVariable::ParetoVariable(double m, double s, double b) : RandomVariable (ParetoVariableImpl (m, s, b)) {} -double -ParetoVariable::GetSingleValue(double m, double s, double b) -{ - return ParetoVariableImpl::GetSingleValue (m, s, b); -} //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -879,14 +687,7 @@ */ virtual double GetValue(); virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param m Mean value for the distribution. - * \param s Shape (alpha) parameter for the distribution. - * \param b Upper limit on returned values - * \return Random number from a distribution specified by m,s, and b - */ - static double GetSingleValue(double m, double s, double b=0); + private: double m_mean; // Mean value of RV double m_alpha; // Shape parameter @@ -906,20 +707,17 @@ double WeibullVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } double exponent = 1.0 / m_alpha; - double r = m_mean * pow( -log(m_generator->RandU01()), exponent); - if (m_bound != 0 && r > m_bound) return m_bound; - return r; + while(1) + { + double r = m_mean * pow( -log(m_generator->RandU01()), exponent); + if (m_bound == 0 || r <= m_bound) return r; + //otherwise, try again + } } RandomVariableBase* WeibullVariableImpl::Copy() const @@ -927,20 +725,6 @@ return new WeibullVariableImpl(*this); } -double WeibullVariableImpl::GetSingleValue(double m, double s, double b/*=0*/) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - double exponent = 1.0 / s; - double r = m * pow( -log(m_static_generator->RandU01()), exponent); - if (b != 0 && r > b) return b; - return r; -} - WeibullVariable::WeibullVariable() : RandomVariable (WeibullVariableImpl ()) {} @@ -953,11 +737,7 @@ WeibullVariable::WeibullVariable(double m, double s, double b) : RandomVariable (WeibullVariableImpl (m, s, b)) {} -double -WeibullVariable::GetSingleValue(double m, double s, double b) -{ - return WeibullVariableImpl::GetSingleValue (m, s, b); -} + //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- // NormalVariableImpl methods @@ -976,7 +756,7 @@ * \brief Construct a normal random variable with specified mean and variance * \param m Mean value * \param v Variance - * \param b Bound. The NormalVariableImpl is bounded within +-bound. + * \param b Bound. The NormalVariableImpl is bounded within +-bound of the mean. */ NormalVariableImpl(double m, double v, double b = INFINITE_VALUE); @@ -987,18 +767,15 @@ */ virtual double GetValue(); virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param m Mean value - * \param v Variance - * \param b Bound. The NormalVariableImpl is bounded within +-bound. - * \return A random number from a distribution specified by m,v, and b. - */ - static double GetSingleValue(double m, double v, double b = INFINITE_VALUE); + + double GetMean (void) const; + double GetVariance (void) const; + double GetBound (void) const; + private: double m_mean; // Mean value of RV double m_variance; // Mean value of RV - double m_bound; // Bound on value (absolute value) + double m_bound; // Bound on value's difference from the mean (absolute value) bool m_nextValid; // True if next valid double m_next; // The algorithm produces two values at a time static bool m_static_nextValid; @@ -1017,20 +794,14 @@ NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c) : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance), - m_bound(c.m_bound), m_nextValid(false) { } + m_bound(c.m_bound) { } double NormalVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } if (m_nextValid) { // use previously generated m_nextValid = false; @@ -1054,15 +825,15 @@ double x1 = m_mean + v1 * y * sqrt(m_variance); //if x1 is in bounds, return it if (fabs(x1-m_mean) <= m_bound) - { - return x1; - } + { + return x1; + } //otherwise try and return m_next if it is valid else if (m_nextValid) - { - m_nextValid = false; - return m_next; - } + { + m_nextValid = false; + return m_next; + } //otherwise, just run this loop again } } @@ -1073,49 +844,22 @@ return new NormalVariableImpl(*this); } -double NormalVariableImpl::GetSingleValue(double m, double v, double b) +double +NormalVariableImpl::GetMean (void) const { - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - if (m_static_nextValid) - { // use previously generated - m_static_nextValid = false; - return m_static_next; - } - while(1) - { // See Simulation Modeling and Analysis p. 466 (Averill Law) - // for algorithm; basically a Box-Muller transform: - // http://en.wikipedia.org/wiki/Box-Muller_transform - double u1 = m_static_generator->RandU01(); - double u2 = m_static_generator->RandU01();; - double v1 = 2 * u1 - 1; - double v2 = 2 * u2 - 1; - double w = v1 * v1 + v2 * v2; - if (w <= 1.0) - { // Got good pair - double y = sqrt((-2 * log(w))/w); - m_static_next = m + v2 * y * sqrt(v); - //if next is in bounds, it is valid - m_static_nextValid = fabs(m_static_next-m) <= b;; - double x1 = m + v1 * y * sqrt(v); - //if x1 is in bounds, return it - if (fabs(x1-m) <= b) - { - return x1; - } - //otherwise try and return m_next if it is valid - else if (m_static_nextValid) - { - m_static_nextValid = false; - return m_static_next; - } - //otherwise, just run this loop again - } - } + return m_mean; +} + +double +NormalVariableImpl::GetVariance (void) const +{ + return m_variance; +} + +double +NormalVariableImpl::GetBound (void) const +{ + return m_bound; } NormalVariable::NormalVariable() @@ -1127,17 +871,6 @@ NormalVariable::NormalVariable(double m, double v, double b) : RandomVariable (NormalVariableImpl (m, v, b)) {} -double -NormalVariable::GetSingleValue(double m, double v) -{ - return NormalVariableImpl::GetSingleValue (m, v); -} -double -NormalVariable::GetSingleValue(double m, double v, double b) -{ - return NormalVariableImpl::GetSingleValue (m, v, b); -} - //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -1200,21 +933,27 @@ double EmpiricalVariableImpl::GetValue() { // Return a value from the empirical distribution // This code based (loosely) on code by Bruce Mah (Thanks Bruce!) - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } - if (emp.size() == 0) return 0.0; // HuH? No empirical data - if (!validated) Validate(); // Insure in non-decreasing + { + m_generator = new RngStream(); + } + if (emp.size() == 0) + { + return 0.0; // HuH? No empirical data + } + if (!validated) + { + Validate(); // Insure in non-decreasing + } double r = m_generator->RandU01(); - if (r <= emp.front().cdf)return emp.front().value; // Less than first - if (r >= emp.back().cdf) return emp.back().value; // Greater than last + if (r <= emp.front().cdf) + { + return emp.front().value; // Less than first + } + if (r >= emp.back().cdf) + { + return emp.back().value; // Greater than last + } // Binary search std::vector::size_type bottom = 0; std::vector::size_type top = emp.size() - 1; @@ -1228,8 +967,14 @@ r); } // Not here, adjust bounds - if (r < emp[c].cdf) top = c - 1; - else bottom = c + 1; + if (r < emp[c].cdf) + { + top = c - 1; + } + else + { + bottom = c + 1; + } } } @@ -1366,7 +1111,10 @@ double DeterministicVariableImpl::GetValue() { - if (next == count) next = 0; + if (next == count) + { + next = 0; + } return data[next++]; } @@ -1395,13 +1143,7 @@ */ virtual double GetValue (); virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param mu mu parameter of the underlying normal distribution - * \param sigma sigma parameter of the underlying normal distribution - * \return A random number from the distribution specified by mu and sigma - */ - static double GetSingleValue(double mu, double sigma); + private: double m_mu; double m_sigma; @@ -1447,16 +1189,10 @@ double LogNormalVariableImpl::GetValue () { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } double u, v, r2, normal, z; do @@ -1478,42 +1214,9 @@ return z; } -double LogNormalVariableImpl::GetSingleValue (double mu, double sigma) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - double u, v, r2, normal, z; - do - { - /* choose x,y in uniform square (-1,-1) to (+1,+1) */ - u = -1 + 2 * m_static_generator->RandU01 (); - v = -1 + 2 * m_static_generator->RandU01 (); - - /* see if it is in the unit circle */ - r2 = u * u + v * v; - } - while (r2 > 1.0 || r2 == 0); - - normal = u * sqrt (-2.0 * log (r2) / r2); - - z = exp (sigma * normal + mu); - - return z; -} - LogNormalVariable::LogNormalVariable (double mu, double sigma) : RandomVariable (LogNormalVariableImpl (mu, sigma)) {} -double -LogNormalVariable::GetSingleValue(double mu, double sigma) -{ - return LogNormalVariableImpl::GetSingleValue (mu, sigma); -} - //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -1542,14 +1245,7 @@ */ virtual double GetValue(); virtual RandomVariableBase* Copy(void) const; -public: - /** - * \param s Low end of the range - * \param l High end of the range - * \param mean mean of the distribution - * \return A triangularly distributed random number between s and l - */ - static double GetSingleValue(double s, double l, double mean); + private: double m_min; double m_max; @@ -1568,21 +1264,19 @@ double TriangularVariableImpl::GetValue() { - if(!RandomVariableBase::initialized) - { - RandomVariableBase::Initialize(); - } if(!m_generator) - { - m_generator = new RngStream(); - m_generator->InitializeStream(); - m_generator->ResetNthSubstream(RandomVariableBase::runNumber); - } + { + m_generator = new RngStream(); + } double u = m_generator->RandU01(); if(u <= (m_mode - m_min) / (m_max - m_min) ) - return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) ); + { + return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) ); + } else - return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) ); + { + return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) ); + } } RandomVariableBase* TriangularVariableImpl::Copy() const @@ -1590,34 +1284,92 @@ return new TriangularVariableImpl(*this); } -double TriangularVariableImpl::GetSingleValue(double s, double l, double mean) -{ - if(!RandomVariableBase::m_static_generator) - { - RandomVariableBase::Initialize(); // sets the static package seed - RandomVariableBase::m_static_generator = new RngStream(); - RandomVariableBase::m_static_generator->InitializeStream(); - } - double mode = 3.0*mean-s-l; - double u = m_static_generator->RandU01(); - if(u <= (mode - s) / (l - s) ) - return s + sqrt(u * (l - s) * (mode - s) ); - else - return l - sqrt( (1-u) * (l - s) * (l - mode) ); -} - TriangularVariable::TriangularVariable() : RandomVariable (TriangularVariableImpl ()) {} TriangularVariable::TriangularVariable(double s, double l, double mean) : RandomVariable (TriangularVariableImpl (s,l,mean)) {} -double -TriangularVariable::GetSingleValue(double s, double l, double mean) + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// ZipfVariableImpl +class ZipfVariableImpl : public RandomVariableBase { +public: + /** + * \param n the number of possible items + * \param alpha the alpha parameter + */ + ZipfVariableImpl (long n, double alpha); + + /** + * \A zipf variable with N=1 and alpha=0 + */ + ZipfVariableImpl (); + + /** + * \return A random value from this distribution + */ + virtual double GetValue (); + virtual RandomVariableBase* Copy(void) const; + +private: + long m_n; + double m_alpha; + double m_c; //the normalization constant +}; + + +RandomVariableBase* ZipfVariableImpl::Copy () const { - return TriangularVariableImpl::GetSingleValue (s,l,mean); + return new ZipfVariableImpl (m_n, m_alpha); } +ZipfVariableImpl::ZipfVariableImpl () + :m_n(1), m_alpha(0), m_c(1) +{ +} + + +ZipfVariableImpl::ZipfVariableImpl (long n, double alpha) + :m_n(n), m_alpha(alpha), m_c(0) +{ + //calculate the normalization constant c + for(int i=1;i<=n;i++) + m_c+=(1.0/pow((double)i,alpha)); + m_c=1.0/m_c; +} + +double +ZipfVariableImpl::GetValue () +{ + if(!m_generator) + { + m_generator = new RngStream(); + } + + double u = m_generator->RandU01(); + double sum_prob=0,zipf_value=0; + for(int i=1;i<=m_n;i++) + { + sum_prob+=m_c/pow((double)i,m_alpha); + if(sum_prob>u) + { + zipf_value=i; + break; + } + } + return zipf_value; +} + +ZipfVariable::ZipfVariable () + : RandomVariable (ZipfVariableImpl ()) +{} + +ZipfVariable::ZipfVariable (long n, double alpha) + : RandomVariable (ZipfVariableImpl (n, alpha)) +{} + std::ostream &operator << (std::ostream &os, const RandomVariable &var) { @@ -1634,6 +1386,17 @@ os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax (); return os; } + NormalVariableImpl *normal = dynamic_cast (base); + if (normal != 0) + { + os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance (); + double bound = normal->GetBound (); + if (bound != NormalVariableImpl::INFINITE_VALUE) + { + os << ":" << bound; + } + return os; + } // XXX: support other distributions os.setstate (std::ios_base::badbit); return os; @@ -1679,6 +1442,44 @@ var = UniformVariable (a, b); } } + else if (type == "Normal") + { + if (value.size () == 0) + { + var = NormalVariable (); + } + else + { + tmp = value.find (":"); + if (tmp == value.npos) + { + NS_FATAL_ERROR ("bad Normal value: " << value); + } + std::string::size_type tmp2; + std::string sub = value.substr (tmp + 1, value.npos); + tmp2 = sub.find (":"); + if (tmp2 == value.npos) + { + istringstream issA (value.substr (0, tmp)); + istringstream issB (sub); + double a, b; + issA >> a; + issB >> b; + var = NormalVariable (a, b); + } + else + { + istringstream issA (value.substr (0, tmp)); + istringstream issB (sub.substr (0, tmp2)); + istringstream issC (sub.substr (tmp2 + 1, value.npos)); + double a, b, c; + issA >> a; + issB >> b; + issC >> c; + var = NormalVariable (a, b, c); + } + } + } else { NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << type); @@ -1691,6 +1492,7 @@ }//namespace ns3 + #ifdef RUN_SELF_TESTS #include "test.h" @@ -1747,47 +1549,26 @@ } } - // Test GetSingleValue - { - vector samples; - const int NSAMPLES = 10000; - double sum = 0; - for (int n = NSAMPLES; n; --n) - { - double value = LogNormalVariable::GetSingleValue (mu, sigma); - sum += value; - samples.push_back (value); - } - double obtained_mean = sum / NSAMPLES; - sum = 0; - for (vector::iterator iter = samples.begin (); iter != samples.end (); iter++) - { - double tmp = (*iter - obtained_mean); - sum += tmp*tmp; - } - double obtained_stddev = sqrt (sum / (NSAMPLES - 1)); - - if (not (obtained_mean/desired_mean > 0.90 and obtained_mean/desired_mean < 1.10)) - { - result = false; - Failure () << "Obtained LogNormalVariable::GetSingleValue mean value " << obtained_mean - << ", expected " << desired_mean << std::endl; - } - - if (not (obtained_stddev/desired_stddev > 0.90 and obtained_stddev/desired_stddev < 1.10)) - { - result = false; - Failure () << "Obtained LogNormalVariable::GetSingleValue stddev value " << obtained_stddev << - ", expected " << desired_stddev << std::endl; - } - } - // Test attribute serialization { - RandomVariableValue val; - val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ()); - RandomVariable rng = val.Get (); - NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2"); + { + RandomVariableValue val; + val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ()); + RandomVariable rng = val.Get (); + NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2"); + } + { + RandomVariableValue val; + val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ()); + RandomVariable rng = val.Get (); + NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2"); + } + { + RandomVariableValue val; + val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker ()); + RandomVariable rng = val.Get (); + NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2:0.15"); + } } return result; --- a/src/core/random-variable.h Wed Feb 04 17:47:29 2009 -0800 +++ a/src/core/random-variable.h Wed May 27 15:48:57 2009 +0200 @@ -16,10 +16,11 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // Author: Rajib Bhattacharjea +// Author: Hadi Arbabi // -#ifndef __random_variable_h__ -#define __random_variable_h__ +#ifndef __random_variable_h +#define __random_variable_h #include #include @@ -39,6 +40,70 @@ class RandomVariableBase; +class SeedManager +{ +public: + + /** + * \brief set the seed + * it will duplicate the seed value 6 times + * \code + * SeedManger::SetSeed(15); + * UniformVariable x(2,3); //these will give the same output everytime + * ExponentialVariable y(120); //as long as the seed stays the same + * \endcode + * \param seed + * + * Note, while the underlying RNG takes six integer values as a seed; + * it is sufficient to set these all to the same integer, so we provide + * a simpler interface here that just takes one integer. + */ + static void SetSeed (uint32_t seed); + + /** + * \brief Get the seed value + * \return the seed value + * + * Note: returns the first of the six seed values used in the underlying RNG + */ + static uint32_t GetSeed (); + + /** + * \brief Set the run number of simulation + * + * \code + * SeedManager::SetSeed(12); + * int N = atol(argv[1]); //read in run number from command line + * SeedManager::SetRun(N); + * UniformVariable x(0,10); + * ExponentialVariable y(2902); + * \endcode + * In this example, N could successivly be equal to 1,2,3, etc. and the user + * would continue to get independent runs out of the single simulation. For + * this simple example, the following might work: + * \code + * ./simulation 0 + * ...Results for run 0:... + * + * ./simulation 1 + * ...Results for run 1:... + * \endcode + */ + static void SetRun (uint32_t run); + /** + * \returns the current run number + * @sa SetRun + */ + static uint32_t GetRun (void); + + /** + * \brief Check if seed value is valid if wanted to be used as seed + * \return true if valid and false if invalid + */ + static bool CheckSeed (uint32_t seed); +}; + + /** * \brief The basic RNG for NS-3. * \ingroup randomvariable @@ -73,95 +138,6 @@ * \return Integer cast of ::GetValue() */ uint32_t GetInteger (void) const; - - /** - * \brief Get the internal state of the RNG - * - * This function is for power users who understand the inner workings - * of the underlying RngStream method used. It returns the internal - * state of the RNG via the input parameter. - * \param seed Output parameter; gets overwritten with the internal state of - * of the RNG. - */ - void GetSeed(uint32_t seed[6]) const; - - /** - * \brief Set seeding behavior - * - * Specify whether the POSIX device /dev/random is to - * be used for seeding. When this is used, the underlying - * generator is seeded with data from /dev/random instead of - * being seeded based upon the time of day. For this to be effective, - * it must be called before the creation of the first instance of a - * RandomVariable or subclass. Example: - * \code - * RandomVariable::UseDevRandom(); - * UniformVariable x(2,3); //these are seeded randomly - * ExponentialVariable y(120); //etc - * \endcode - * \param udr True if /dev/random desired. - */ - static void UseDevRandom(bool udr = true); - - /** - * \brief Use the global seed to force precisely reproducible results. - * - * It is often desirable to create a simulation that uses random - * numbers, while at the same time is completely reproducible. - * Specifying this set of six random seeds initializes the - * random number generator with the specified seed. - * Once this is set, all generators will produce fixed output - * from run to run. This is because each time a new generator is created, - * the underlying RngStream deterministically creates a new seed based upon - * the old one, hence a "stream" of RNGs. Example: - * \code - * RandomVariable::UseGlobalSeed(...); - * UniformVariable x(2,3); //these will give the same output everytime - * ExponentialVariable y(120); //as long as the seed stays the same - * \endcode - * \param s0 - * \param s1 - * \param s2 - * \param s3 - * \param s4 - * \param s5 - * \return True if seed is valid. - */ - static void UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, - uint32_t s3, uint32_t s4, uint32_t s5); - - /** - * \brief Set the run number of this simulation - * - * These RNGs have the ability to give independent sets of trials for a fixed - * global seed. For example, suppose one sets up a simulation with - * RandomVariables with a given global seed. Suppose the user wanted to - * retry the same simulation with different random values for validity, - * statistical rigor, etc. The user could either change the global seed and - * re-run the simulation, or could use this facility to increment all of the - * RNGs to a next substream state. This predictably advances the internal - * state of all RandomVariables n steps. This should be called immediately - * after the global seed is set, and before the creation of any - * RandomVariables. For example: - * \code - * RandomVariable::UseGlobalSeed(1,2,3,4,5,6); - * int N = atol(argv[1]); //read in run number from command line - * RandomVariable::SetRunNumber(N); - * UniformVariable x(0,10); - * ExponentialVariable y(2902); - * \endcode - * In this example, N could successivly be equal to 1,2,3, etc. and the user - * would continue to get independent runs out of the single simulation. For - * this simple example, the following might work: - * \code - * ./simulation 0 - * ...Results for run 0:... - * - * ./simulation 1 - * ...Results for run 1:... - * \endcode - */ - static void SetRunNumber(uint32_t n); private: friend std::ostream &operator << (std::ostream &os, const RandomVariable &var); @@ -204,13 +180,33 @@ * \param l High end of the range */ UniformVariable(double s, double l); -public: + /** + * \brief call RandomVariable::GetValue + * \return A floating point random value + * + * Note: we have to re-implement this method here because the method is + * overloaded below for the two-argument variant and the c++ name resolution + * rules don't work well with overloads split between parent and child + * classes. + */ + double GetValue (void) const; + + /** + * \brief Returns a random double with the specified range + * \param s Low end of the range + * \param l High end of the range + * \return A floating point random value + */ + double GetValue(double s, double l); + + /** + * \brief Returns a random unsigned integer from the interval [s,l] including both ends. * \param s Low end of the range * \param l High end of the range - * \return A uniformly distributed random number between s and l + * \return A random unsigned integer value. */ - static double GetSingleValue(double s, double l); + uint32_t GetInteger (uint32_t s, uint32_t l); }; /** @@ -334,12 +330,6 @@ */ ExponentialVariable(double m, double b); - /** - * \param m The mean of the distribution from which the return value is drawn - * \param b The upper bound value desired, beyond which values get clipped - * \return A random number from an exponential distribution with mean m - */ - static double GetSingleValue(double m, double b=0); }; /** @@ -402,16 +392,6 @@ */ ParetoVariable(double m, double s, double b); - /** - * \param m The mean value of the distribution from which the return value - * is drawn. - * \param s The shape parameter of the distribution from which the return - * value is drawn. - * \param b The upper bound to which to restrict return values - * \return A random number from a Pareto distribution with mean m and shape - * parameter s. - */ - static double GetSingleValue(double m, double s, double b=0); }; /** @@ -464,13 +444,7 @@ * \param b Upper limit on returned values */ WeibullVariable(double m, double s, double b); - /** - * \param m Mean value for the distribution. - * \param s Shape (alpha) parameter for the distribution. - * \param b Upper limit on returned values - * \return Random number from a distribution specified by m,s, and b - */ - static double GetSingleValue(double m, double s, double b=0); + }; /** @@ -511,22 +485,6 @@ * [mean-bound,mean+bound] */ NormalVariable(double m, double v, double b); - - /** - * \param m Mean value - * \param v Variance - * \return A random number from a distribution specified by m, and v. - */ - static double GetSingleValue(double m, double v); - - /** - * \param m Mean value - * \param v Variance - * \param b Bound. The NormalVariable is bounded symetrically about the mean - * [mean-bound,mean+bound] - * \return A random number from a distribution specified by m,v, and b. - */ - static double GetSingleValue(double m, double v, double b); }; /** @@ -633,13 +591,25 @@ * \param sigma sigma parameter of the lognormal distribution */ LogNormalVariable (double mu, double sigma); +}; +/** + * \brief Zipf Distributed random var (between 1 and n included) + * \ingroup randomvariable + * + */ +class ZipfVariable : public RandomVariable +{ +public: /** - * \param mu mu parameter of the underlying normal distribution - * \param sigma sigma parameter of the underlying normal distribution - * \return A random number from the distribution specified by mu and sigma + * \param n the number of possible items + * \param alpha the alpha parameter */ - static double GetSingleValue(double mu, double sigma); + ZipfVariable (long n, double alpha); + /** + * A zipf variable with N=1 and alpha=0 + */ + ZipfVariable (); }; /** @@ -666,13 +636,7 @@ * \param mean mean of the distribution */ TriangularVariable(double s, double l, double mean); - /** - * \param s Low end of the range - * \param l High end of the range - * \param mean mean of the distribution - * \return A triangularly distributed random number between s and l - */ - static double GetSingleValue(double s, double l, double mean); + }; std::ostream &operator << (std::ostream &os, const RandomVariable &var);