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