src/core/random-variable.cc
author Mathieu Lacage <mathieu.lacage@sophia.inria.fr>
Sat, 04 Jul 2009 08:15:48 +0200
changeset 4654 2eaebe77d66b
parent 4317 0a250f44e0ed
permissions -rw-r--r--
Added tag ns-3.5 for changeset c975274c9707
     1 /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
     2 //
     3 // Copyright (c) 2006 Georgia Tech Research Corporation
     4 //
     5 // This program is free software; you can redistribute it and/or modify
     6 // it under the terms of the GNU General Public License version 2 as
     7 // published by the Free Software Foundation;
     8 //
     9 // This program is distributed in the hope that it will be useful,
    10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
    11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    12 // GNU General Public License for more details.
    13 //
    14 // You should have received a copy of the GNU General Public License
    15 // along with this program; if not, write to the Free Software
    16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
    17 //
    18 // Author: Rajib Bhattacharjea<raj.b@gatech.edu>
    19 // Author: Hadi Arbabi<marbabi@cs.odu.edu>
    20 //
    21 
    22 #include <iostream>
    23 
    24 #include <math.h>
    25 #include <stdlib.h>
    26 #include <sys/time.h>			// for gettimeofday
    27 #include <unistd.h>
    28 #include <iostream>
    29 #include <sys/types.h>
    30 #include <sys/stat.h>
    31 #include <fcntl.h>       
    32 #include <sstream>
    33 
    34 
    35 #include "assert.h"
    36 #include "config.h"
    37 #include "integer.h"
    38 #include "random-variable.h"
    39 #include "rng-stream.h"
    40 #include "fatal-error.h"
    41 
    42 using namespace std;
    43 
    44 namespace ns3{
    45 
    46 //-----------------------------------------------------------------------------
    47 // Seed Manager
    48 //-----------------------------------------------------------------------------
    49 
    50 uint32_t SeedManager::GetSeed()
    51 {
    52   uint32_t s[6];
    53   RngStream::GetPackageSeed (s);
    54   NS_ASSERT(
    55               s[0] == s[1] &&
    56               s[0] == s[2] &&
    57               s[0] == s[3] &&
    58               s[0] == s[4] &&
    59               s[0] == s[5]    
    60             );
    61   return s[0];
    62 }
    63 
    64 void SeedManager::SetSeed(uint32_t seed)
    65 {
    66   Config::SetGlobal("RngSeed", IntegerValue(seed));
    67 }
    68 
    69 void SeedManager::SetRun(uint32_t run)
    70 {
    71   Config::SetGlobal("RngRun", IntegerValue(run));
    72 }
    73 
    74 uint32_t SeedManager::GetRun()
    75 {
    76   return RngStream::GetPackageRun ();
    77 }
    78 
    79 bool SeedManager::CheckSeed (uint32_t seed)
    80 {
    81   return RngStream::CheckSeed(seed);
    82 }
    83 
    84 //-----------------------------------------------------------------------------
    85 //-----------------------------------------------------------------------------
    86 // RandomVariableBase methods
    87 
    88 
    89 class RandomVariableBase 
    90 {
    91 public:
    92   RandomVariableBase ();
    93   RandomVariableBase (const RandomVariableBase &o);
    94   virtual ~RandomVariableBase();
    95   virtual double  GetValue() = 0;
    96   virtual uint32_t GetInteger();
    97   virtual RandomVariableBase*   Copy(void) const = 0;
    98 
    99 protected:
   100   RngStream* m_generator;  //underlying generator being wrapped
   101 };
   102 
   103 RandomVariableBase::RandomVariableBase() 
   104   : m_generator(NULL)
   105 {
   106 }
   107 
   108 RandomVariableBase::RandomVariableBase(const RandomVariableBase& r)
   109   :m_generator(0)
   110 {
   111   if (r.m_generator)
   112     {
   113       m_generator = new RngStream(*r.m_generator);
   114     }
   115 }
   116 
   117 RandomVariableBase::~RandomVariableBase()
   118 {
   119   delete m_generator;
   120 }
   121 
   122 uint32_t RandomVariableBase::GetInteger() 
   123 {
   124   return (uint32_t)GetValue();
   125 }
   126 
   127 //-------------------------------------------------------
   128 
   129 RandomVariable::RandomVariable()
   130   : m_variable (0)
   131 {}
   132 RandomVariable::RandomVariable(const RandomVariable&o)
   133   : m_variable (o.m_variable->Copy ())
   134 {}
   135 RandomVariable::RandomVariable (const RandomVariableBase &variable)
   136   : m_variable (variable.Copy ())
   137 {}
   138 RandomVariable &
   139 RandomVariable::operator = (const RandomVariable &o)
   140 {
   141   if (&o == this)
   142     {
   143       return *this;
   144     }
   145   delete m_variable;
   146   m_variable = o.m_variable->Copy ();
   147   return *this;
   148 }
   149 RandomVariable::~RandomVariable()
   150 {
   151   delete m_variable;
   152 }
   153 double  
   154 RandomVariable::GetValue (void) const
   155 {
   156   return m_variable->GetValue ();
   157 }
   158 
   159 uint32_t 
   160 RandomVariable::GetInteger (void) const
   161 {
   162   return m_variable->GetInteger ();
   163 }
   164 
   165 RandomVariableBase *
   166 RandomVariable::Peek (void) const
   167 {
   168   return m_variable;
   169 }
   170 
   171 
   172 ATTRIBUTE_VALUE_IMPLEMENT (RandomVariable);
   173 ATTRIBUTE_CHECKER_IMPLEMENT (RandomVariable);
   174 
   175 //-----------------------------------------------------------------------------
   176 //-----------------------------------------------------------------------------
   177 // UniformVariableImpl
   178 
   179 class UniformVariableImpl : public RandomVariableBase {
   180 public:
   181   /**
   182    * Creates a uniform random number generator in the
   183    * range [0.0 .. 1.0).
   184    */
   185   UniformVariableImpl();
   186 
   187   /**
   188    * Creates a uniform random number generator with the specified range
   189    * \param s Low end of the range
   190    * \param l High end of the range
   191    */
   192   UniformVariableImpl(double s, double l);
   193 
   194   UniformVariableImpl(const UniformVariableImpl& c);
   195 
   196   double GetMin (void) const;
   197   double GetMax (void) const;
   198   
   199   /**
   200    * \return A value between low and high values specified by the constructor
   201    */
   202   virtual double GetValue();
   203 
   204   /**
   205    * \return A value between low and high values specified by parameters
   206    */
   207   virtual double GetValue(double s, double l);
   208   
   209   virtual RandomVariableBase*  Copy(void) const;
   210 
   211 private:
   212   double m_min;
   213   double m_max;
   214 };
   215 
   216 UniformVariableImpl::UniformVariableImpl() 
   217   : m_min(0), m_max(1.0) { }
   218   
   219 UniformVariableImpl::UniformVariableImpl(double s, double l) 
   220   : m_min(s), m_max(l) { }
   221 
   222 UniformVariableImpl::UniformVariableImpl(const UniformVariableImpl& c) 
   223   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max) { }
   224 
   225 double 
   226 UniformVariableImpl::GetMin (void) const
   227 {
   228   return m_min;
   229 }
   230 double 
   231 UniformVariableImpl::GetMax (void) const
   232 {
   233   return m_max;
   234 }
   235 
   236 
   237 double UniformVariableImpl::GetValue()
   238 {
   239   if(!m_generator)
   240     {
   241       m_generator = new RngStream();
   242     }
   243   return m_min + m_generator->RandU01() * (m_max - m_min);
   244 }
   245 
   246 double UniformVariableImpl::GetValue(double s, double l) 
   247 {
   248   if(!m_generator)
   249     {
   250       m_generator = new RngStream();
   251     }
   252     return s + m_generator->RandU01() * (l-s); 
   253 }
   254 
   255 RandomVariableBase* UniformVariableImpl::Copy() const
   256 {
   257   return new UniformVariableImpl(*this);
   258 }
   259 
   260 UniformVariable::UniformVariable()
   261   : RandomVariable (UniformVariableImpl ())
   262 {}
   263 UniformVariable::UniformVariable(double s, double l)
   264   : RandomVariable (UniformVariableImpl (s, l))
   265 {}
   266 
   267 double UniformVariable::GetValue(void) const
   268 {
   269   return this->RandomVariable::GetValue ();
   270 }
   271 
   272 double UniformVariable::GetValue(double s, double l)
   273 {
   274   return ((UniformVariableImpl*)Peek())->GetValue(s,l);
   275 }
   276 
   277 uint32_t UniformVariable::GetInteger (uint32_t s, uint32_t l)
   278 {
   279   NS_ASSERT(s <= l);
   280   return static_cast<uint32_t>( GetValue(s, l+1) );
   281 }
   282 
   283 //-----------------------------------------------------------------------------
   284 //-----------------------------------------------------------------------------
   285 // ConstantVariableImpl methods
   286 
   287 class ConstantVariableImpl : public RandomVariableBase { 
   288 
   289 public:
   290   /**
   291    * Construct a ConstantVariableImpl RNG that returns zero every sample
   292    */
   293   ConstantVariableImpl();
   294   
   295   /**
   296    * Construct a ConstantVariableImpl RNG that returns the specified value
   297    * every sample.
   298    * \param c Unchanging value for this RNG.
   299    */
   300   ConstantVariableImpl(double c);
   301 
   302 
   303   ConstantVariableImpl(const ConstantVariableImpl& c) ;
   304 
   305   /**
   306    * \brief Specify a new constant RNG for this generator.
   307    * \param c New constant value for this RNG.
   308    */
   309   void    NewConstant(double c);
   310 
   311   /**
   312    * \return The constant value specified
   313    */
   314   virtual double  GetValue();
   315   virtual uint32_t GetInteger();
   316   virtual RandomVariableBase*   Copy(void) const;
   317 private:
   318   double m_const;
   319 };
   320 
   321 ConstantVariableImpl::ConstantVariableImpl() 
   322   : m_const(0) { }
   323 
   324 ConstantVariableImpl::ConstantVariableImpl(double c) 
   325   : m_const(c) { };
   326   
   327 ConstantVariableImpl::ConstantVariableImpl(const ConstantVariableImpl& c) 
   328   : RandomVariableBase(c), m_const(c.m_const) { }
   329 
   330 void ConstantVariableImpl::NewConstant(double c) 
   331   { m_const = c;}
   332   
   333 double ConstantVariableImpl::GetValue()
   334 {
   335   return m_const;
   336 }
   337 
   338 uint32_t ConstantVariableImpl::GetInteger()
   339 {
   340   return (uint32_t)m_const;
   341 }
   342 
   343 RandomVariableBase* ConstantVariableImpl::Copy() const
   344 {
   345   return new ConstantVariableImpl(*this);
   346 }
   347 
   348 ConstantVariable::ConstantVariable()
   349   : RandomVariable (ConstantVariableImpl ())
   350 {}
   351 ConstantVariable::ConstantVariable(double c)
   352   : RandomVariable (ConstantVariableImpl (c))
   353 {}
   354 void 
   355 ConstantVariable::SetConstant(double c)
   356 {
   357   *this = ConstantVariable (c);
   358 }
   359 
   360 //-----------------------------------------------------------------------------
   361 //-----------------------------------------------------------------------------
   362 // SequentialVariableImpl methods
   363 
   364 
   365 class SequentialVariableImpl : public RandomVariableBase {
   366 
   367 public:
   368   /**
   369    * \brief Constructor for the SequentialVariableImpl RNG.
   370    *
   371    * The four parameters define the sequence.  For example
   372    * SequentialVariableImpl(0,5,1,2) creates a RNG that has the sequence
   373    * 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0 ...
   374    * \param f First value of the sequence.
   375    * \param l One more than the last value of the sequence.
   376    * \param i Increment between sequence values
   377    * \param c Number of times each member of the sequence is repeated
   378    */
   379   SequentialVariableImpl(double f, double l, double i = 1, uint32_t c = 1);
   380 
   381   /**
   382    * \brief Constructor for the SequentialVariableImpl RNG.
   383    *
   384    * Differs from the first only in that the increment parameter is a
   385    * random variable
   386    * \param f First value of the sequence.
   387    * \param l One more than the last value of the sequence.
   388    * \param i Reference to a RandomVariableBase for the sequence increment
   389    * \param c Number of times each member of the sequence is repeated
   390    */
   391   SequentialVariableImpl(double f, double l, const RandomVariable& i, uint32_t c = 1);
   392 
   393   SequentialVariableImpl(const SequentialVariableImpl& c);
   394   
   395   ~SequentialVariableImpl();
   396   /**
   397    * \return The next value in the Sequence
   398    */
   399   virtual double GetValue();
   400   virtual RandomVariableBase*  Copy(void) const;
   401 private:
   402   double m_min;
   403   double m_max;
   404   RandomVariable  m_increment;
   405   uint32_t  m_consecutive;
   406   double m_current;
   407   uint32_t  m_currentConsecutive;
   408 };
   409 
   410 SequentialVariableImpl::SequentialVariableImpl(double f, double l, double i, uint32_t c)
   411   : m_min(f), m_max(l), m_increment(ConstantVariable(i)), m_consecutive(c),
   412     m_current(f), m_currentConsecutive(0)
   413 {}
   414 
   415 SequentialVariableImpl::SequentialVariableImpl(double f, double l, const RandomVariable& i, uint32_t c)
   416   : m_min(f), m_max(l), m_increment(i), m_consecutive(c),
   417     m_current(f), m_currentConsecutive(0)
   418 {}
   419 
   420 SequentialVariableImpl::SequentialVariableImpl(const SequentialVariableImpl& c)
   421   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max),
   422     m_increment(c.m_increment), m_consecutive(c.m_consecutive),
   423     m_current(c.m_current), m_currentConsecutive(c.m_currentConsecutive)
   424 {}
   425 
   426 SequentialVariableImpl::~SequentialVariableImpl()
   427 {}
   428 
   429 double SequentialVariableImpl::GetValue()
   430 { // Return a sequential series of values
   431   double r = m_current;
   432   if (++m_currentConsecutive == m_consecutive)
   433     { // Time to advance to next
   434       m_currentConsecutive = 0;
   435       m_current += m_increment.GetValue();
   436       if (m_current >= m_max)
   437         m_current = m_min + (m_current - m_max);
   438     }
   439   return r;
   440 }
   441 
   442 RandomVariableBase* SequentialVariableImpl::Copy() const
   443 {
   444   return new SequentialVariableImpl(*this);
   445 }
   446 
   447 SequentialVariable::SequentialVariable(double f, double l, double i, uint32_t c)
   448   : RandomVariable (SequentialVariableImpl (f, l, i, c))
   449 {}
   450 SequentialVariable::SequentialVariable(double f, double l, const RandomVariable& i, uint32_t c)
   451   : RandomVariable (SequentialVariableImpl (f, l, i, c))
   452 {}
   453 
   454 //-----------------------------------------------------------------------------
   455 //-----------------------------------------------------------------------------
   456 // ExponentialVariableImpl methods
   457 
   458 class ExponentialVariableImpl : public RandomVariableBase { 
   459 public:
   460   /**
   461    * Constructs an exponential random variable  with a mean
   462    * value of 1.0.
   463    */
   464   ExponentialVariableImpl();
   465 
   466   /**
   467    * \brief Constructs an exponential random variable with a specified mean
   468    * \param m Mean value for the random variable
   469    */
   470   explicit ExponentialVariableImpl(double m);
   471 
   472   /**
   473    * \brief Constructs an exponential random variable with spefified
   474    * \brief mean and upper limit.
   475    *
   476    * Since exponential distributions can theoretically return unbounded values,
   477    * it is sometimes useful to specify a fixed upper limit.  Note however when
   478    * the upper limit is specified, the true mean of the distribution is 
   479    * slightly smaller than the mean value specified.
   480    * \param m Mean value of the random variable
   481    * \param b Upper bound on returned values
   482    */
   483   ExponentialVariableImpl(double m, double b);
   484 
   485   ExponentialVariableImpl(const ExponentialVariableImpl& c);
   486   
   487   /**
   488    * \return A random value from this exponential distribution
   489    */
   490   virtual double GetValue();
   491   virtual RandomVariableBase* Copy(void) const;
   492 
   493 private:
   494   double m_mean;  // Mean value of RV
   495   double m_bound; // Upper bound on value (if non-zero)
   496 };
   497 
   498 ExponentialVariableImpl::ExponentialVariableImpl() 
   499   : m_mean(1.0), m_bound(0) { }
   500   
   501 ExponentialVariableImpl::ExponentialVariableImpl(double m) 
   502   : m_mean(m), m_bound(0) { }
   503   
   504 ExponentialVariableImpl::ExponentialVariableImpl(double m, double b) 
   505   : m_mean(m), m_bound(b) { }
   506   
   507 ExponentialVariableImpl::ExponentialVariableImpl(const ExponentialVariableImpl& c) 
   508   : RandomVariableBase(c), m_mean(c.m_mean), m_bound(c.m_bound) { }
   509 
   510 double ExponentialVariableImpl::GetValue()
   511 {
   512   if(!m_generator)
   513     {
   514       m_generator = new RngStream();
   515     }
   516   while(1)
   517     {
   518       double r = -m_mean*log(m_generator->RandU01());
   519       if (m_bound == 0 || r <= m_bound) return r;
   520       //otherwise, try again
   521     }
   522 }
   523 
   524 RandomVariableBase* ExponentialVariableImpl::Copy() const
   525 {
   526   return new ExponentialVariableImpl(*this);
   527 }
   528 
   529 ExponentialVariable::ExponentialVariable()
   530   : RandomVariable (ExponentialVariableImpl ())
   531 {}
   532 ExponentialVariable::ExponentialVariable(double m)
   533   : RandomVariable (ExponentialVariableImpl (m))
   534 {}
   535 ExponentialVariable::ExponentialVariable(double m, double b)
   536   : RandomVariable (ExponentialVariableImpl (m, b))
   537 {}
   538 
   539 //-----------------------------------------------------------------------------
   540 //-----------------------------------------------------------------------------
   541 // ParetoVariableImpl methods
   542 class ParetoVariableImpl : public RandomVariableBase {
   543 public:
   544   /**
   545    * Constructs a pareto random variable with a mean of 1 and a shape
   546    * parameter of 1.5
   547    */
   548   ParetoVariableImpl();
   549 
   550   /**
   551    * Constructs a pareto random variable with specified mean and shape
   552    * parameter of 1.5
   553    * \param m Mean value of the distribution
   554    */
   555   explicit ParetoVariableImpl(double m);
   556 
   557   /**
   558    * Constructs a pareto random variable with the specified mean value and
   559    * shape parameter.
   560    * \param m Mean value of the distribution
   561    * \param s Shape parameter for the distribution
   562    */
   563   ParetoVariableImpl(double m, double s);
   564 
   565   /**
   566    * \brief Constructs a pareto random variable with the specified mean
   567    * \brief value, shape (alpha), and upper bound.
   568    *
   569    * Since pareto distributions can theoretically return unbounded values,
   570    * it is sometimes useful to specify a fixed upper limit.  Note however
   571    * when the upper limit is specified, the true mean of the distribution
   572    * is slightly smaller than the mean value specified.
   573    * \param m Mean value
   574    * \param s Shape parameter
   575    * \param b Upper limit on returned values
   576    */
   577   ParetoVariableImpl(double m, double s, double b);
   578 
   579   ParetoVariableImpl(const ParetoVariableImpl& c);
   580   
   581   /**
   582    * \return A random value from this Pareto distribution
   583    */
   584   virtual double GetValue();
   585   virtual RandomVariableBase* Copy() const;
   586 
   587 private:
   588   double m_mean;  // Mean value of RV
   589   double m_shape; // Shape parameter
   590   double m_bound; // Upper bound on value (if non-zero)
   591 };
   592 
   593 ParetoVariableImpl::ParetoVariableImpl() 
   594   : m_mean(1.0), m_shape(1.5), m_bound(0) { }
   595 
   596 ParetoVariableImpl::ParetoVariableImpl(double m) 
   597   : m_mean(m), m_shape(1.5), m_bound(0) { }
   598 
   599 ParetoVariableImpl::ParetoVariableImpl(double m, double s) 
   600     : m_mean(m), m_shape(s), m_bound(0) { }
   601 
   602 ParetoVariableImpl::ParetoVariableImpl(double m, double s, double b) 
   603   : m_mean(m), m_shape(s), m_bound(b) { }
   604 
   605 ParetoVariableImpl::ParetoVariableImpl(const ParetoVariableImpl& c) 
   606   : RandomVariableBase(c), m_mean(c.m_mean), m_shape(c.m_shape), 
   607     m_bound(c.m_bound) { }
   608 
   609 double ParetoVariableImpl::GetValue()
   610 {
   611   if(!m_generator)
   612     {
   613       m_generator = new RngStream();
   614     }
   615   double scale = m_mean * ( m_shape - 1.0) / m_shape;
   616   while(1)
   617     {
   618       double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
   619       if (m_bound == 0 || r <= m_bound) return r;
   620       //otherwise, try again
   621     }
   622 }
   623 
   624 RandomVariableBase* ParetoVariableImpl::Copy() const
   625 {
   626   return new ParetoVariableImpl(*this);
   627 }
   628 
   629 ParetoVariable::ParetoVariable ()
   630   : RandomVariable (ParetoVariableImpl ())
   631 {}
   632 ParetoVariable::ParetoVariable(double m)
   633   : RandomVariable (ParetoVariableImpl (m))
   634 {}
   635 ParetoVariable::ParetoVariable(double m, double s)
   636   : RandomVariable (ParetoVariableImpl (m, s))
   637 {}
   638 ParetoVariable::ParetoVariable(double m, double s, double b)
   639   : RandomVariable (ParetoVariableImpl (m, s, b))
   640 {}
   641 
   642 //-----------------------------------------------------------------------------
   643 //-----------------------------------------------------------------------------
   644 // WeibullVariableImpl methods
   645 
   646 class WeibullVariableImpl : public RandomVariableBase {
   647 public:
   648   /**
   649    * Constructs a weibull random variable  with a mean
   650    * value of 1.0 and a shape (alpha) parameter of 1
   651    */
   652   WeibullVariableImpl();
   653 
   654 
   655   /**
   656    * Constructs a weibull random variable with the specified mean
   657    * value and a shape (alpha) parameter of 1.5.
   658    * \param m mean value of the distribution
   659    */
   660    WeibullVariableImpl(double m) ;
   661 
   662   /**
   663    * Constructs a weibull random variable with the specified mean
   664    * value and a shape (alpha).
   665    * \param m Mean value for the distribution.
   666    * \param s Shape (alpha) parameter for the distribution.
   667    */
   668   WeibullVariableImpl(double m, double s);
   669 
   670    /**
   671    * \brief Constructs a weibull random variable with the specified mean
   672    * \brief value, shape (alpha), and upper bound.
   673    * Since WeibullVariableImpl distributions can theoretically return unbounded values,
   674    * it is sometimes usefull to specify a fixed upper limit.  Note however
   675    * that when the upper limit is specified, the true mean of the distribution
   676    * is slightly smaller than the mean value specified.
   677    * \param m Mean value for the distribution.
   678    * \param s Shape (alpha) parameter for the distribution.
   679    * \param b Upper limit on returned values
   680    */
   681   WeibullVariableImpl(double m, double s, double b);
   682 
   683   WeibullVariableImpl(const WeibullVariableImpl& c);
   684   
   685   /**
   686    * \return A random value from this Weibull distribution
   687    */
   688   virtual double GetValue();
   689   virtual RandomVariableBase* Copy(void) const;
   690 
   691 private:
   692   double m_mean;  // Mean value of RV
   693   double m_alpha; // Shape parameter
   694   double m_bound; // Upper bound on value (if non-zero)
   695 };
   696 
   697 WeibullVariableImpl::WeibullVariableImpl() : m_mean(1.0), m_alpha(1), m_bound(0) { }
   698 WeibullVariableImpl::WeibullVariableImpl(double m) 
   699   : m_mean(m), m_alpha(1), m_bound(0) { }
   700 WeibullVariableImpl::WeibullVariableImpl(double m, double s) 
   701   : m_mean(m), m_alpha(s), m_bound(0) { }
   702 WeibullVariableImpl::WeibullVariableImpl(double m, double s, double b) 
   703   : m_mean(m), m_alpha(s), m_bound(b) { };
   704 WeibullVariableImpl::WeibullVariableImpl(const WeibullVariableImpl& c) 
   705   : RandomVariableBase(c), m_mean(c.m_mean), m_alpha(c.m_alpha),
   706     m_bound(c.m_bound) { }
   707 
   708 double WeibullVariableImpl::GetValue()
   709 {
   710   if(!m_generator)
   711     {
   712       m_generator = new RngStream();
   713     }
   714   double exponent = 1.0 / m_alpha;
   715   while(1)
   716     {
   717       double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
   718       if (m_bound == 0 || r <= m_bound) return r;
   719       //otherwise, try again
   720     }
   721 }
   722 
   723 RandomVariableBase* WeibullVariableImpl::Copy() const
   724 {
   725   return new WeibullVariableImpl(*this);
   726 }
   727 
   728 WeibullVariable::WeibullVariable()
   729   : RandomVariable (WeibullVariableImpl ())
   730 {}
   731 WeibullVariable::WeibullVariable(double m)
   732   : RandomVariable (WeibullVariableImpl (m))
   733 {}
   734 WeibullVariable::WeibullVariable(double m, double s)
   735   : RandomVariable (WeibullVariableImpl (m, s))
   736 {}
   737 WeibullVariable::WeibullVariable(double m, double s, double b)
   738   : RandomVariable (WeibullVariableImpl (m, s, b))
   739 {}
   740 
   741 //-----------------------------------------------------------------------------
   742 //-----------------------------------------------------------------------------
   743 // NormalVariableImpl methods
   744 
   745 class NormalVariableImpl : public RandomVariableBase { // Normally Distributed random var
   746 
   747 public:
   748    static const double INFINITE_VALUE;
   749   /**
   750    * Constructs an normal random variable  with a mean
   751    * value of 0 and variance of 1.
   752    */ 
   753   NormalVariableImpl();
   754 
   755   /**
   756    * \brief Construct a normal random variable with specified mean and variance
   757    * \param m Mean value
   758    * \param v Variance
   759    * \param b Bound.  The NormalVariableImpl is bounded within +-bound of the mean.
   760    */ 
   761   NormalVariableImpl(double m, double v, double b = INFINITE_VALUE);
   762 
   763   NormalVariableImpl(const NormalVariableImpl& c);
   764   
   765   /**
   766    * \return A value from this normal distribution
   767    */
   768   virtual double GetValue();
   769   virtual RandomVariableBase* Copy(void) const;
   770 
   771   double GetMean (void) const;
   772   double GetVariance (void) const;
   773   double GetBound (void) const;
   774 
   775 private:
   776   double m_mean;      // Mean value of RV
   777   double m_variance;  // Mean value of RV
   778   double m_bound;     // Bound on value's difference from the mean (absolute value)
   779   bool   m_nextValid; // True if next valid
   780   double m_next;      // The algorithm produces two values at a time
   781   static bool   m_static_nextValid;
   782   static double m_static_next;
   783 };
   784 
   785 bool         NormalVariableImpl::m_static_nextValid = false;
   786 double       NormalVariableImpl::m_static_next;
   787 const double NormalVariableImpl::INFINITE_VALUE = 1e307;
   788 
   789 NormalVariableImpl::NormalVariableImpl() 
   790   : m_mean(0.0), m_variance(1.0), m_bound(INFINITE_VALUE), m_nextValid(false){}
   791 
   792 NormalVariableImpl::NormalVariableImpl(double m, double v, double b/*=INFINITE_VALUE*/)
   793   : m_mean(m), m_variance(v), m_bound(b), m_nextValid(false) { }
   794 
   795 NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c)
   796   : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance),
   797     m_bound(c.m_bound) { }
   798 
   799 double NormalVariableImpl::GetValue()
   800 {
   801   if(!m_generator)
   802     {
   803       m_generator = new RngStream();
   804     }
   805   if (m_nextValid)
   806     { // use previously generated
   807       m_nextValid = false;
   808       return m_next;
   809     }
   810   while(1)
   811     { // See Simulation Modeling and Analysis p. 466 (Averill Law)
   812       // for algorithm; basically a Box-Muller transform:
   813       // http://en.wikipedia.org/wiki/Box-Muller_transform
   814       double u1 = m_generator->RandU01();
   815       double u2 = m_generator->RandU01();;
   816       double v1 = 2 * u1 - 1;
   817       double v2 = 2 * u2 - 1;
   818       double w = v1 * v1 + v2 * v2;
   819       if (w <= 1.0)
   820         { // Got good pair
   821           double y = sqrt((-2 * log(w))/w);
   822           m_next = m_mean + v2 * y * sqrt(m_variance);
   823           //if next is in bounds, it is valid
   824           m_nextValid = fabs(m_next-m_mean) <= m_bound;
   825           double x1 = m_mean + v1 * y * sqrt(m_variance);
   826           //if x1 is in bounds, return it
   827           if (fabs(x1-m_mean) <= m_bound)
   828             {
   829               return x1;
   830             }
   831           //otherwise try and return m_next if it is valid
   832           else if (m_nextValid)
   833 	    {
   834 	      m_nextValid = false;
   835 	      return m_next;
   836 	    }
   837           //otherwise, just run this loop again
   838         }
   839     }
   840 }
   841 
   842 RandomVariableBase* NormalVariableImpl::Copy() const
   843 {
   844   return new NormalVariableImpl(*this);
   845 }
   846 
   847 double
   848 NormalVariableImpl::GetMean (void) const
   849 {
   850   return m_mean;
   851 }
   852 
   853 double
   854 NormalVariableImpl::GetVariance (void) const
   855 {
   856   return m_variance;
   857 }
   858 
   859 double
   860 NormalVariableImpl::GetBound (void) const
   861 {
   862   return m_bound;
   863 }
   864 
   865 NormalVariable::NormalVariable()
   866   : RandomVariable (NormalVariableImpl ())
   867 {}
   868 NormalVariable::NormalVariable(double m, double v)
   869   : RandomVariable (NormalVariableImpl (m, v))
   870 {}
   871 NormalVariable::NormalVariable(double m, double v, double b)
   872   : RandomVariable (NormalVariableImpl (m, v, b))
   873 {}
   874 
   875 //-----------------------------------------------------------------------------
   876 //-----------------------------------------------------------------------------
   877 class EmpiricalVariableImpl : public RandomVariableBase {
   878 public:
   879   /**
   880    * Constructor for the EmpiricalVariableImpl random variables.
   881    */
   882   explicit EmpiricalVariableImpl();
   883 
   884   virtual ~EmpiricalVariableImpl();
   885   EmpiricalVariableImpl(const EmpiricalVariableImpl& c);
   886   /**
   887    * \return A value from this empirical distribution
   888    */
   889   virtual double GetValue();
   890   virtual RandomVariableBase* Copy(void) const;
   891   /**
   892    * \brief Specifies a point in the empirical distribution
   893    * \param v The function value for this point
   894    * \param c Probability that the function is less than or equal to v
   895    */
   896   virtual void CDF(double v, double c);  // Value, prob <= Value
   897 
   898 private:
   899   class ValueCDF {
   900   public:
   901     ValueCDF();
   902     ValueCDF(double v, double c);
   903     ValueCDF(const ValueCDF& c);
   904     double value;
   905     double    cdf;
   906   };
   907   virtual void Validate();  // Insure non-decreasing emiprical values
   908   virtual double Interpolate(double, double, double, double, double);
   909   bool validated; // True if non-decreasing validated
   910   std::vector<ValueCDF> emp;       // Empicical CDF
   911 };
   912 
   913 
   914 // ValueCDF methods
   915 EmpiricalVariableImpl::ValueCDF::ValueCDF() 
   916   : value(0.0), cdf(0.0){ }
   917 EmpiricalVariableImpl::ValueCDF::ValueCDF(double v, double c) 
   918   : value(v), cdf(c) { }
   919 EmpiricalVariableImpl::ValueCDF::ValueCDF(const ValueCDF& c) 
   920   : value(c.value), cdf(c.cdf) { }
   921 
   922 //-----------------------------------------------------------------------------
   923 //-----------------------------------------------------------------------------
   924 // EmpiricalVariableImpl methods
   925 EmpiricalVariableImpl::EmpiricalVariableImpl() 
   926   : validated(false) { }
   927 
   928 EmpiricalVariableImpl::EmpiricalVariableImpl(const EmpiricalVariableImpl& c)
   929   : RandomVariableBase(c), validated(c.validated), emp(c.emp) { }
   930 
   931 EmpiricalVariableImpl::~EmpiricalVariableImpl() { }
   932 
   933 double EmpiricalVariableImpl::GetValue()
   934 { // Return a value from the empirical distribution
   935   // This code based (loosely) on code by Bruce Mah (Thanks Bruce!)
   936   if(!m_generator)
   937     {
   938       m_generator = new RngStream();
   939     }
   940   if (emp.size() == 0) 
   941     {
   942       return 0.0; // HuH? No empirical data
   943     }
   944   if (!validated) 
   945     {
   946       Validate();      // Insure in non-decreasing
   947     }
   948   double r = m_generator->RandU01();
   949   if (r <= emp.front().cdf)
   950     {
   951       return emp.front().value; // Less than first
   952     }
   953   if (r >= emp.back().cdf) 
   954     { 
   955       return emp.back().value;  // Greater than last
   956     }
   957   // Binary search
   958   std::vector<ValueCDF>::size_type bottom = 0;
   959   std::vector<ValueCDF>::size_type top = emp.size() - 1;
   960   while(1)
   961     {
   962       std::vector<ValueCDF>::size_type c = (top + bottom) / 2;
   963       if (r >= emp[c].cdf && r < emp[c+1].cdf)
   964         { // Found it
   965           return Interpolate(emp[c].cdf, emp[c+1].cdf,
   966                              emp[c].value, emp[c+1].value,
   967                              r);
   968         }
   969       // Not here, adjust bounds
   970       if (r < emp[c].cdf)
   971         {
   972           top    = c - 1;
   973         }
   974       else
   975         {
   976           bottom = c + 1;
   977         }
   978     }
   979 }
   980 
   981 RandomVariableBase* EmpiricalVariableImpl::Copy() const
   982 {
   983   return new EmpiricalVariableImpl(*this);
   984 }
   985 
   986 void EmpiricalVariableImpl::CDF(double v, double c)
   987 { // Add a new empirical datapoint to the empirical cdf
   988   // NOTE.   These MUST be inserted in non-decreasing order
   989   emp.push_back(ValueCDF(v, c));
   990 }
   991 
   992 void EmpiricalVariableImpl::Validate()
   993 {
   994   ValueCDF prior;
   995   for (std::vector<ValueCDF>::size_type i = 0; i < emp.size(); ++i)
   996     {
   997       ValueCDF& current = emp[i];
   998       if (current.value < prior.value || current.cdf < prior.cdf)
   999         { // Error
  1000           cerr << "Empirical Dist error,"
  1001                << " current value " << current.value
  1002                << " prior value "   << prior.value
  1003                << " current cdf "   << current.cdf
  1004                << " prior cdf "     << prior.cdf << endl;
  1005           NS_FATAL_ERROR("Empirical Dist error");
  1006         }
  1007       prior = current;
  1008     }
  1009   validated = true;
  1010 }
  1011 
  1012 double EmpiricalVariableImpl::Interpolate(double c1, double c2,
  1013                                 double v1, double v2, double r)
  1014 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)
  1015   return (v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
  1016 }
  1017 
  1018 EmpiricalVariable::EmpiricalVariable()
  1019   : RandomVariable (EmpiricalVariableImpl ())
  1020 {}
  1021 EmpiricalVariable::EmpiricalVariable (const RandomVariableBase &variable)
  1022   : RandomVariable (variable)
  1023 {}
  1024 void 
  1025 EmpiricalVariable::CDF(double v, double c)
  1026 {
  1027   EmpiricalVariableImpl *impl = dynamic_cast<EmpiricalVariableImpl *> (Peek ());
  1028   NS_ASSERT (impl);
  1029   impl->CDF (v, c);
  1030 }
  1031 
  1032 
  1033 //-----------------------------------------------------------------------------
  1034 //-----------------------------------------------------------------------------
  1035 // IntegerValue EmpiricalVariableImpl methods
  1036 class IntEmpiricalVariableImpl : public EmpiricalVariableImpl {
  1037 public:
  1038 
  1039   IntEmpiricalVariableImpl();
  1040   
  1041   virtual RandomVariableBase* Copy(void) const;
  1042   /**
  1043    * \return An integer value from this empirical distribution
  1044    */
  1045   virtual uint32_t GetInteger();
  1046 private:
  1047   virtual double Interpolate(double, double, double, double, double);
  1048 };
  1049 
  1050 
  1051 IntEmpiricalVariableImpl::IntEmpiricalVariableImpl() { }
  1052 
  1053 uint32_t IntEmpiricalVariableImpl::GetInteger()
  1054 {
  1055   return (uint32_t)GetValue();
  1056 }
  1057 
  1058 RandomVariableBase* IntEmpiricalVariableImpl::Copy() const
  1059 {
  1060   return new IntEmpiricalVariableImpl(*this);
  1061 }
  1062 
  1063 double IntEmpiricalVariableImpl::Interpolate(double c1, double c2,
  1064                                    double v1, double v2, double r)
  1065 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)
  1066   return ceil(v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
  1067 }
  1068 
  1069 IntEmpiricalVariable::IntEmpiricalVariable()
  1070   : EmpiricalVariable (IntEmpiricalVariableImpl ())
  1071 {}
  1072 
  1073 //-----------------------------------------------------------------------------
  1074 //-----------------------------------------------------------------------------
  1075 // DeterministicVariableImpl
  1076 class DeterministicVariableImpl : public RandomVariableBase 
  1077 {
  1078 
  1079 public:
  1080   /**
  1081    * \brief Constructor
  1082    *
  1083    * Creates a generator that returns successive elements of the d array
  1084    * on successive calls to ::Value().  Note that the d pointer is copied
  1085    * for use by the generator (shallow-copy), not its contents, so the 
  1086    * contents of the array d points to have to remain unchanged for the use 
  1087    * of DeterministicVariableImpl to be meaningful.
  1088    * \param d Pointer to array of random values to return in sequence
  1089    * \param c Number of values in the array
  1090    */
  1091   explicit DeterministicVariableImpl(double* d, uint32_t c);
  1092 
  1093   virtual ~DeterministicVariableImpl();
  1094   /**
  1095    * \return The next value in the deterministic sequence
  1096    */
  1097   virtual double GetValue();
  1098   virtual RandomVariableBase* Copy(void) const;
  1099 private:
  1100   uint32_t   count;
  1101   uint32_t   next;
  1102   double* data;
  1103 };
  1104 
  1105 DeterministicVariableImpl::DeterministicVariableImpl(double* d, uint32_t c)
  1106     : count(c), next(c), data(d)
  1107 { // Nothing else needed
  1108 }
  1109 
  1110 DeterministicVariableImpl::~DeterministicVariableImpl() { }
  1111   
  1112 double DeterministicVariableImpl::GetValue()
  1113 {
  1114   if (next == count) 
  1115     {
  1116       next = 0;
  1117     }
  1118   return data[next++];
  1119 }
  1120 
  1121 RandomVariableBase* DeterministicVariableImpl::Copy() const
  1122 {
  1123   return new DeterministicVariableImpl(*this);
  1124 }
  1125 
  1126 DeterministicVariable::DeterministicVariable(double* d, uint32_t c)
  1127   : RandomVariable (DeterministicVariableImpl (d, c))
  1128 {}
  1129 
  1130 //-----------------------------------------------------------------------------
  1131 //-----------------------------------------------------------------------------
  1132 // LogNormalVariableImpl
  1133 class LogNormalVariableImpl : public RandomVariableBase { 
  1134 public:
  1135   /**
  1136    * \param mu mu parameter of the lognormal distribution
  1137    * \param sigma sigma parameter of the lognormal distribution
  1138    */
  1139   LogNormalVariableImpl (double mu, double sigma);
  1140 
  1141   /**
  1142    * \return A random value from this distribution
  1143    */
  1144   virtual double GetValue ();
  1145   virtual RandomVariableBase* Copy(void) const;
  1146 
  1147 private:
  1148   double m_mu;
  1149   double m_sigma;
  1150 };
  1151 
  1152 
  1153 RandomVariableBase* LogNormalVariableImpl::Copy () const
  1154 {
  1155   return new LogNormalVariableImpl (*this);
  1156 }
  1157 
  1158 LogNormalVariableImpl::LogNormalVariableImpl (double mu, double sigma)
  1159     :m_mu(mu), m_sigma(sigma) 
  1160 {
  1161 }
  1162 
  1163 // The code from this function was adapted from the GNU Scientific
  1164 // Library 1.8:
  1165 /* randist/lognormal.c
  1166  * 
  1167  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  1168  * 
  1169  * This program is free software; you can redistribute it and/or modify
  1170  * it under the terms of the GNU General Public License as published by
  1171  * the Free Software Foundation; either version 2 of the License, or (at
  1172  * your option) any later version.
  1173  * 
  1174  * This program is distributed in the hope that it will be useful, but
  1175  * WITHOUT ANY WARRANTY; without even the implied warranty of
  1176  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  1177  * General Public License for more details.
  1178  * 
  1179  * You should have received a copy of the GNU General Public License
  1180  * along with this program; if not, write to the Free Software
  1181  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  1182  */
  1183 /* The lognormal distribution has the form 
  1184 
  1185    p(x) dx = 1/(x * sqrt(2 pi sigma^2)) exp(-(ln(x) - zeta)^2/2 sigma^2) dx
  1186 
  1187    for x > 0. Lognormal random numbers are the exponentials of
  1188    gaussian random numbers */
  1189 double
  1190 LogNormalVariableImpl::GetValue ()
  1191 {
  1192   if(!m_generator)
  1193     {
  1194       m_generator = new RngStream();
  1195     }
  1196   double u, v, r2, normal, z;
  1197 
  1198   do
  1199     {
  1200       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
  1201 
  1202       u = -1 + 2 * m_generator->RandU01 ();
  1203       v = -1 + 2 * m_generator->RandU01 ();
  1204 
  1205       /* see if it is in the unit circle */
  1206       r2 = u * u + v * v;
  1207     }
  1208   while (r2 > 1.0 || r2 == 0);
  1209 
  1210   normal = u * sqrt (-2.0 * log (r2) / r2);
  1211 
  1212   z =  exp (m_sigma * normal + m_mu);
  1213 
  1214   return z;
  1215 }
  1216 
  1217 LogNormalVariable::LogNormalVariable (double mu, double sigma)
  1218   : RandomVariable (LogNormalVariableImpl (mu, sigma))
  1219 {}
  1220 
  1221 //-----------------------------------------------------------------------------
  1222 //-----------------------------------------------------------------------------
  1223 // GammaVariableImpl
  1224 class GammaVariableImpl : public RandomVariableBase
  1225 {
  1226 public:
  1227   /**
  1228    * \param alpha alpha parameter of the gamma distribution
  1229    * \param beta beta parameter of the gamma distribution
  1230    */
  1231   GammaVariableImpl (double alpha, double beta);
  1232 
  1233   /**
  1234    * \return A random value from this distribution
  1235    */
  1236   virtual double GetValue ();
  1237 
  1238   /**
  1239    * \return A random value from the gamma distribution with parameters alpha
  1240    * and beta
  1241    */
  1242   double GetValue(double alpha, double beta);
  1243 
  1244   virtual RandomVariableBase* Copy(void) const;
  1245 
  1246 private:
  1247   double m_alpha;
  1248   double m_beta;
  1249   NormalVariable m_normal;
  1250 };
  1251 
  1252 
  1253 RandomVariableBase* GammaVariableImpl::Copy () const
  1254 {
  1255   return new GammaVariableImpl (m_alpha, m_beta);
  1256 }
  1257 
  1258 GammaVariableImpl::GammaVariableImpl (double alpha, double beta)
  1259   : m_alpha(alpha), m_beta(beta) 
  1260 {
  1261 }
  1262 
  1263 double
  1264 GammaVariableImpl::GetValue ()
  1265 {
  1266   return GetValue(m_alpha, m_beta);
  1267 }
  1268 
  1269 /*
  1270   The code for the following generator functions was adapted from ns-2
  1271   tools/ranvar.cc
  1272 
  1273   Originally the algorithm was devised by Marsaglia in 2000:
  1274   G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
  1275   ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
  1276 
  1277   The Gamma distribution density function has the form 
  1278 
  1279 	                     x^(alpha-1) * exp(-x/beta)
  1280 	p(x; alpha, beta) = ----------------------------
  1281                              beta^alpha * Gamma(alpha)
  1282 
  1283   for x > 0.
  1284 */
  1285 double
  1286 GammaVariableImpl::GetValue (double alpha, double beta)
  1287 {
  1288   if(!m_generator)
  1289     {
  1290       m_generator = new RngStream();
  1291     }
  1292 
  1293   if (alpha < 1)
  1294     {
  1295       double u = m_generator->RandU01 ();
  1296       return GetValue(1.0 + alpha, beta) * pow (u, 1.0 / alpha);
  1297     }
  1298 	
  1299   double x, v, u;
  1300   double d = alpha - 1.0 / 3.0;
  1301   double c = (1.0 / 3.0) / sqrt (d);
  1302 
  1303   while (1)
  1304     {
  1305       do
  1306         {
  1307           x = m_normal.GetValue ();
  1308           v = 1.0 + c * x;
  1309         } while (v <= 0);
  1310 
  1311       v = v * v * v;
  1312       u = m_generator->RandU01 ();
  1313       if (u < 1 - 0.0331 * x * x * x * x)
  1314         {
  1315           break;
  1316         }
  1317       if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
  1318         {
  1319           break;
  1320         }
  1321     }
  1322 
  1323   return beta * d * v;
  1324 }
  1325 
  1326 GammaVariable::GammaVariable ()
  1327   : RandomVariable (GammaVariableImpl (1.0, 1.0))
  1328 {
  1329 }
  1330 
  1331 GammaVariable::GammaVariable (double alpha, double beta)
  1332   : RandomVariable (GammaVariableImpl (alpha, beta))
  1333 {
  1334 }
  1335 
  1336 double GammaVariable::GetValue(void) const
  1337 {
  1338   return this->RandomVariable::GetValue ();
  1339 }
  1340 
  1341 double GammaVariable::GetValue(double alpha, double beta) const
  1342 {
  1343   return ((GammaVariableImpl*)Peek())->GetValue(alpha, beta);
  1344 }
  1345 
  1346 //-----------------------------------------------------------------------------
  1347 //-----------------------------------------------------------------------------
  1348 // ErlangVariableImpl
  1349 
  1350 class ErlangVariableImpl : public RandomVariableBase
  1351 {
  1352 public:
  1353   /**
  1354    * \param k k parameter of the Erlang distribution
  1355    * \param lambda lambda parameter of the Erlang distribution
  1356    */
  1357   ErlangVariableImpl (unsigned int k, double lambda);
  1358 
  1359   /**
  1360    * \return A random value from this distribution
  1361    */
  1362   virtual double GetValue ();
  1363 
  1364   /**
  1365    * \return A random value from the Erlang distribution with parameters k and
  1366    * lambda.
  1367    */
  1368   double GetValue(unsigned int k, double lambda);
  1369 
  1370   virtual RandomVariableBase* Copy(void) const;
  1371 
  1372 private:
  1373   unsigned int m_k;
  1374   double m_lambda;
  1375 };
  1376 
  1377 
  1378 RandomVariableBase* ErlangVariableImpl::Copy () const
  1379 {
  1380   return new ErlangVariableImpl (m_k, m_lambda);
  1381 }
  1382 
  1383 ErlangVariableImpl::ErlangVariableImpl (unsigned int k, double lambda)
  1384   : m_k(k), m_lambda(lambda)
  1385 {
  1386 }
  1387 
  1388 double
  1389 ErlangVariableImpl::GetValue ()
  1390 {
  1391   return GetValue(m_k, m_lambda);
  1392 }
  1393 
  1394 /*
  1395   The code for the following generator functions was adapted from ns-2
  1396   tools/ranvar.cc
  1397 
  1398   The Erlang distribution density function has the form 
  1399 
  1400 	                   x^(k-1) * exp(-x/lambda)
  1401 	p(x; k, lambda) = ---------------------------
  1402                              lambda^k * (k-1)!
  1403 
  1404   for x > 0.
  1405 */
  1406 double
  1407 ErlangVariableImpl::GetValue (unsigned int k, double lambda)
  1408 {
  1409   if(!m_generator)
  1410     {
  1411       m_generator = new RngStream();
  1412     }
  1413 
  1414   ExponentialVariable exponential(lambda);
  1415 
  1416   double result = 0;
  1417   for (unsigned int i = 0; i < k; ++i)
  1418     {
  1419       result += exponential.GetValue();
  1420     }
  1421 
  1422   return result;
  1423 }
  1424 
  1425 ErlangVariable::ErlangVariable ()
  1426   : RandomVariable (ErlangVariableImpl (1, 1.0))
  1427 {
  1428 }
  1429 
  1430 ErlangVariable::ErlangVariable (unsigned int k, double lambda)
  1431   : RandomVariable (ErlangVariableImpl (k, lambda))
  1432 {
  1433 }
  1434 
  1435 double ErlangVariable::GetValue(void) const
  1436 {
  1437   return this->RandomVariable::GetValue ();
  1438 }
  1439 
  1440 double ErlangVariable::GetValue(unsigned int k, double lambda) const
  1441 {
  1442   return ((ErlangVariableImpl*)Peek())->GetValue(k, lambda);
  1443 }
  1444 
  1445 //-----------------------------------------------------------------------------
  1446 //-----------------------------------------------------------------------------
  1447 // TriangularVariableImpl methods
  1448 class TriangularVariableImpl : public RandomVariableBase {
  1449 public:
  1450   /**
  1451    * Creates a triangle distribution random number generator in the
  1452    * range [0.0 .. 1.0), with mean of 0.5
  1453    */
  1454   TriangularVariableImpl();
  1455 
  1456   /**
  1457    * Creates a triangle distribution random number generator with the specified
  1458    * range
  1459    * \param s Low end of the range
  1460    * \param l High end of the range
  1461    * \param mean mean of the distribution
  1462    */
  1463   TriangularVariableImpl(double s, double l, double mean);
  1464 
  1465   TriangularVariableImpl(const TriangularVariableImpl& c);
  1466   
  1467   /**
  1468    * \return A value from this distribution
  1469    */
  1470   virtual double GetValue();
  1471   virtual RandomVariableBase*  Copy(void) const;
  1472 
  1473 private:
  1474   double m_min;
  1475   double m_max;
  1476   double m_mode;  //easier to work with the mode internally instead of the mean
  1477                   //they are related by the simple: mean = (min+max+mode)/3
  1478 };
  1479 
  1480 TriangularVariableImpl::TriangularVariableImpl() 
  1481   : m_min(0), m_max(1), m_mode(0.5) { }
  1482   
  1483 TriangularVariableImpl::TriangularVariableImpl(double s, double l, double mean) 
  1484   : m_min(s), m_max(l), m_mode(3.0*mean-s-l) { }
  1485   
  1486 TriangularVariableImpl::TriangularVariableImpl(const TriangularVariableImpl& c) 
  1487   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max), m_mode(c.m_mode) { }
  1488 
  1489 double TriangularVariableImpl::GetValue()
  1490 {
  1491   if(!m_generator)
  1492     {
  1493       m_generator = new RngStream();
  1494     }
  1495   double u = m_generator->RandU01();
  1496   if(u <= (m_mode - m_min) / (m_max - m_min) )
  1497     {
  1498       return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) );
  1499     }
  1500   else
  1501     {
  1502       return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) );
  1503     }
  1504 }
  1505 
  1506 RandomVariableBase* TriangularVariableImpl::Copy() const
  1507 {
  1508   return new TriangularVariableImpl(*this);
  1509 }
  1510 
  1511 TriangularVariable::TriangularVariable()
  1512   : RandomVariable (TriangularVariableImpl ())
  1513 {}
  1514 TriangularVariable::TriangularVariable(double s, double l, double mean)
  1515   : RandomVariable (TriangularVariableImpl (s,l,mean))
  1516 {}
  1517 
  1518 //-----------------------------------------------------------------------------
  1519 //-----------------------------------------------------------------------------
  1520 // ZipfVariableImpl
  1521 class ZipfVariableImpl : public RandomVariableBase { 
  1522 public:
  1523   /**
  1524    * \param n the number of possible items
  1525    * \param alpha the alpha parameter
  1526    */
  1527   ZipfVariableImpl (long n, double alpha);
  1528 
  1529   /**
  1530    * \A zipf variable with N=1 and alpha=0
  1531    */
  1532   ZipfVariableImpl ();
  1533 
  1534   /**
  1535    * \return A random value from this distribution
  1536    */
  1537   virtual double GetValue ();
  1538   virtual RandomVariableBase* Copy(void) const;
  1539 
  1540 private:
  1541   long m_n;
  1542   double m_alpha;
  1543   double m_c; //the normalization constant
  1544 };
  1545 
  1546 
  1547 RandomVariableBase* ZipfVariableImpl::Copy () const
  1548 {
  1549   return new ZipfVariableImpl (m_n, m_alpha);
  1550 }
  1551 
  1552 ZipfVariableImpl::ZipfVariableImpl ()
  1553     :m_n(1), m_alpha(0), m_c(1)
  1554 {
  1555 }
  1556 
  1557 
  1558 ZipfVariableImpl::ZipfVariableImpl (long n, double alpha)
  1559     :m_n(n), m_alpha(alpha), m_c(0)
  1560 {
  1561   //calculate the normalization constant c
  1562   for(int i=1;i<=n;i++)
  1563     {
  1564       m_c+=(1.0/pow((double)i,alpha));
  1565     }
  1566   m_c=1.0/m_c;
  1567 }
  1568 
  1569 double
  1570 ZipfVariableImpl::GetValue ()
  1571 {
  1572   if(!m_generator)
  1573     {
  1574       m_generator = new RngStream();
  1575     }
  1576 
  1577   double u = m_generator->RandU01();
  1578   double sum_prob=0,zipf_value=0;
  1579   for(int i=1;i<=m_n;i++)
  1580     {
  1581       sum_prob+=m_c/pow((double)i,m_alpha);
  1582       if(sum_prob>u)
  1583         {
  1584           zipf_value=i;
  1585           break;
  1586         }
  1587     }
  1588   return zipf_value;
  1589 }
  1590 
  1591 ZipfVariable::ZipfVariable ()
  1592   : RandomVariable (ZipfVariableImpl ())
  1593 {}
  1594 
  1595 ZipfVariable::ZipfVariable (long n, double alpha)
  1596   : RandomVariable (ZipfVariableImpl (n, alpha))
  1597 {}
  1598 
  1599 
  1600 std::ostream &operator << (std::ostream &os, const RandomVariable &var)
  1601 {
  1602   RandomVariableBase *base = var.Peek ();
  1603   ConstantVariableImpl *constant = dynamic_cast<ConstantVariableImpl *> (base);
  1604   if (constant != 0)
  1605     {
  1606       os << "Constant:" << constant->GetValue ();
  1607       return os;
  1608     }
  1609   UniformVariableImpl *uniform = dynamic_cast<UniformVariableImpl *> (base);
  1610   if (uniform != 0)
  1611     {
  1612       os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax ();
  1613       return os;
  1614     }
  1615   NormalVariableImpl *normal = dynamic_cast<NormalVariableImpl *> (base);
  1616   if (normal != 0)
  1617     {
  1618       os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance ();
  1619       double bound = normal->GetBound ();
  1620       if (bound != NormalVariableImpl::INFINITE_VALUE)
  1621         {
  1622           os << ":" << bound;
  1623         }
  1624       return os;
  1625     }
  1626   // XXX: support other distributions
  1627   os.setstate (std::ios_base::badbit);
  1628   return os;
  1629 }
  1630 std::istream &operator >> (std::istream &is, RandomVariable &var)
  1631 {
  1632   std::string value;
  1633   is >> value;
  1634   std::string::size_type tmp;
  1635   tmp = value.find (":");
  1636   if (tmp == std::string::npos)
  1637     {
  1638       is.setstate (std::ios_base::badbit);
  1639       return is;
  1640     }
  1641   std::string type = value.substr (0, tmp);
  1642   value = value.substr (tmp + 1, value.npos);
  1643   if (type == "Constant")
  1644     {
  1645       istringstream iss (value);
  1646       double constant;
  1647       iss >> constant;
  1648       var = ConstantVariable (constant);
  1649     }
  1650   else if (type == "Uniform")
  1651     {
  1652       if (value.size () == 0)
  1653         {
  1654           var = UniformVariable ();
  1655         }
  1656       else
  1657         {
  1658           tmp = value.find (":");
  1659           if (tmp == value.npos)
  1660             {
  1661               NS_FATAL_ERROR ("bad Uniform value: " << value);
  1662             }
  1663           istringstream issA (value.substr (0, tmp));
  1664           istringstream issB (value.substr (tmp + 1, value.npos));
  1665           double a, b;
  1666           issA >> a;
  1667           issB >> b;
  1668           var = UniformVariable (a, b);
  1669         }
  1670     }
  1671   else if (type == "Normal")
  1672     {
  1673       if (value.size () == 0)
  1674         {
  1675           var = NormalVariable ();
  1676         }
  1677       else
  1678         {
  1679           tmp = value.find (":");
  1680           if (tmp == value.npos)
  1681             {
  1682               NS_FATAL_ERROR ("bad Normal value: " << value);
  1683             }
  1684           std::string::size_type tmp2;
  1685           std::string sub = value.substr (tmp + 1, value.npos);
  1686           tmp2 = sub.find (":");
  1687           if (tmp2 == value.npos)
  1688             {
  1689               istringstream issA (value.substr (0, tmp));
  1690               istringstream issB (sub);
  1691               double a, b;
  1692               issA >> a;
  1693               issB >> b;
  1694               var = NormalVariable (a, b);
  1695             }
  1696           else
  1697             {
  1698               istringstream issA (value.substr (0, tmp));
  1699               istringstream issB (sub.substr (0, tmp2));
  1700               istringstream issC (sub.substr (tmp2 + 1, value.npos));
  1701               double a, b, c;
  1702               issA >> a;
  1703               issB >> b;
  1704               issC >> c;
  1705               var = NormalVariable (a, b, c);
  1706             }
  1707         }
  1708     }
  1709   else
  1710     {
  1711       NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << type);
  1712       // XXX: support other distributions.
  1713     }
  1714   return is;
  1715 }
  1716 
  1717 
  1718 
  1719 }//namespace ns3
  1720 
  1721  
  1722 
  1723 #ifdef RUN_SELF_TESTS
  1724 #include "test.h"
  1725 #include <vector>
  1726 
  1727 namespace ns3 {
  1728 
  1729 
  1730 class RandomVariableTest : public Test
  1731 {
  1732 public:
  1733   RandomVariableTest () : Test ("RandomVariable") {}
  1734   virtual bool RunTests (void)
  1735   {
  1736     bool result = true;
  1737     const double desired_mean = 1.0;
  1738     const double desired_stddev = 1.0;
  1739     double tmp = log (1 + (desired_stddev/desired_mean)*(desired_stddev/desired_mean));
  1740     double sigma = sqrt (tmp);
  1741     double mu = log (desired_mean) - 0.5*tmp;
  1742 
  1743     // Test a custom lognormal instance
  1744     {
  1745       LogNormalVariable lognormal (mu, sigma);
  1746       vector<double> samples;
  1747       const int NSAMPLES = 10000;
  1748       double sum = 0;
  1749       for (int n = NSAMPLES; n; --n)
  1750         {
  1751           double value = lognormal.GetValue ();
  1752           sum += value;
  1753           samples.push_back (value);
  1754         }
  1755       double obtained_mean = sum / NSAMPLES;
  1756       sum = 0;
  1757       for (vector<double>::iterator iter = samples.begin (); iter != samples.end (); iter++)
  1758         {
  1759           double tmp = (*iter - obtained_mean);
  1760           sum += tmp*tmp;
  1761         }
  1762       double obtained_stddev = sqrt (sum / (NSAMPLES - 1));
  1763 
  1764       if (not (obtained_mean/desired_mean > 0.90 and obtained_mean/desired_mean < 1.10))
  1765         {
  1766           result = false;
  1767           Failure () << "Obtained lognormal mean value " << obtained_mean << ", expected " << desired_mean << std::endl;
  1768         }
  1769 
  1770       if (not (obtained_stddev/desired_stddev > 0.90 and obtained_stddev/desired_stddev < 1.10))
  1771         {
  1772           result = false;
  1773           Failure () << "Obtained lognormal stddev value " << obtained_stddev <<
  1774             ", expected " << desired_stddev << std::endl;
  1775         }
  1776     }
  1777 
  1778     // Test attribute serialization
  1779     {
  1780       {
  1781         RandomVariableValue val;
  1782         val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ());
  1783         RandomVariable rng = val.Get ();
  1784         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2");
  1785       }
  1786       {
  1787         RandomVariableValue val;
  1788         val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ());
  1789         RandomVariable rng = val.Get ();
  1790         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2");
  1791       }
  1792       {
  1793         RandomVariableValue val;
  1794         val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker ());
  1795         RandomVariable rng = val.Get ();
  1796         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2:0.15");
  1797       }
  1798     }
  1799 
  1800     return result;
  1801   }
  1802 };
  1803 
  1804 
  1805 static RandomVariableTest g_random_variable_tests;
  1806 
  1807 }//namespace ns3
  1808 
  1809 #endif /* RUN_SELF_TESTS */