tomh@346: /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
tomh@346: //
tomh@346: // Copyright (c) 2006 Georgia Tech Research Corporation
tomh@346: //
tomh@346: // This program is free software; you can redistribute it and/or modify
tomh@346: // it under the terms of the GNU General Public License version 2 as
tomh@346: // published by the Free Software Foundation;
tomh@346: //
tomh@346: // This program is distributed in the hope that it will be useful,
tomh@346: // but WITHOUT ANY WARRANTY; without even the implied warranty of
tomh@346: // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
tomh@346: // GNU General Public License for more details.
tomh@346: //
tomh@346: // You should have received a copy of the GNU General Public License
tomh@346: // along with this program; if not, write to the Free Software
tomh@346: // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
tomh@346: //
tomh@346: // Author: Rajib Bhattacharjea<raj.b@gatech.edu>
raj@4218: // Author: Hadi Arbabi<marbabi@cs.odu.edu>
tomh@346: //
tomh@346: 
tomh@346: #include <iostream>
tomh@346: 
tomh@346: #include <math.h>
tomh@346: #include <stdlib.h>
tomh@346: #include <sys/time.h>			// for gettimeofday
tomh@346: #include <unistd.h>
tomh@346: #include <iostream>
tomh@346: #include <sys/types.h>
tomh@346: #include <sys/stat.h>
tomh@346: #include <fcntl.h>       
mathieu@2384: #include <sstream>
tomh@346: 
tomh@346: 
mathieu@2336: #include "assert.h"
raj@4235: #include "config.h"
raj@4234: #include "integer.h"
tomh@346: #include "random-variable.h"
tomh@346: #include "rng-stream.h"
tomh@346: #include "fatal-error.h"
tomh@346: 
tomh@346: using namespace std;
tomh@346: 
tomh@346: namespace ns3{
tomh@346: 
tomh@346: //-----------------------------------------------------------------------------
raj@4218: // Seed Manager
raj@4218: //-----------------------------------------------------------------------------
raj@4218: 
raj@4234: uint32_t SeedManager::GetSeed()
raj@4218: {
raj@4234:   uint32_t s[6];
raj@4234:   RngStream::GetPackageSeed (s);
raj@4234:   NS_ASSERT(
raj@4234:               s[0] == s[1] &&
raj@4234:               s[0] == s[2] &&
raj@4234:               s[0] == s[3] &&
raj@4234:               s[0] == s[4] &&
raj@4234:               s[0] == s[5]    
raj@4234:             );
raj@4234:   return s[0];
raj@4218: }
raj@4218: 
raj@4218: void SeedManager::SetSeed(uint32_t seed)
raj@4218: {
raj@4235:   Config::SetGlobal("RngSeed", IntegerValue(seed));
raj@4218: }
raj@4218: 
raj@4218: void SeedManager::SetRun(uint32_t run)
raj@4218: {
raj@4235:   Config::SetGlobal("RngRun", IntegerValue(run));
raj@4218: }
raj@4218: 
raj@4218: uint32_t SeedManager::GetRun()
raj@4218: {
raj@4224:   return RngStream::GetPackageRun ();
raj@4218: }
raj@4218: 
raj@4218: bool SeedManager::CheckSeed (uint32_t seed)
raj@4218: {
raj@4224:   return RngStream::CheckSeed(seed);
raj@4218: }
raj@4218: 
raj@4218: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // RandomVariableBase methods
mathieu@2336: 
mathieu@2336: 
mathieu@2336: class RandomVariableBase 
mathieu@2336: {
mathieu@2336: public:
mathieu@2336:   RandomVariableBase ();
mathieu@2336:   RandomVariableBase (const RandomVariableBase &o);
mathieu@2336:   virtual ~RandomVariableBase();
mathieu@2336:   virtual double  GetValue() = 0;
mathieu@2439:   virtual uint32_t GetInteger();
mathieu@2336:   virtual RandomVariableBase*   Copy(void) const = 0;
tomh@346: 
mathieu@2336: protected:
mathieu@2336:   RngStream* m_generator;  //underlying generator being wrapped
mathieu@2336: };
mathieu@2336: 
mathieu@2336: RandomVariableBase::RandomVariableBase() 
raj@1806:   : m_generator(NULL)
tomh@346: {
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase::RandomVariableBase(const RandomVariableBase& r)
raj@1811:   :m_generator(0)
raj@360: {
raj@4224:   if (r.m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream(*r.m_generator);
raj@4224:     }
raj@360: }
raj@360: 
mathieu@2336: RandomVariableBase::~RandomVariableBase()
tomh@346: {
tomh@346:   delete m_generator;
tomh@346: }
tomh@346: 
mathieu@2439: uint32_t RandomVariableBase::GetInteger() 
tomh@346: {
tomh@346:   return (uint32_t)GetValue();
tomh@346: }
tomh@346: 
raj@4218: //-------------------------------------------------------
mathieu@2336: 
mathieu@2336: RandomVariable::RandomVariable()
mathieu@2336:   : m_variable (0)
mathieu@2336: {}
mathieu@2336: RandomVariable::RandomVariable(const RandomVariable&o)
mathieu@2336:   : m_variable (o.m_variable->Copy ())
mathieu@2336: {}
mathieu@2336: RandomVariable::RandomVariable (const RandomVariableBase &variable)
mathieu@2336:   : m_variable (variable.Copy ())
mathieu@2336: {}
mathieu@2336: RandomVariable &
mathieu@2336: RandomVariable::operator = (const RandomVariable &o)
mathieu@2336: {
mathieu@2384:   if (&o == this)
mathieu@2384:     {
mathieu@2384:       return *this;
mathieu@2384:     }
mathieu@2336:   delete m_variable;
mathieu@2336:   m_variable = o.m_variable->Copy ();
mathieu@2336:   return *this;
mathieu@2336: }
mathieu@2336: RandomVariable::~RandomVariable()
mathieu@2336: {
mathieu@2336:   delete m_variable;
mathieu@2336: }
mathieu@2336: double  
mathieu@2336: RandomVariable::GetValue (void) const
mathieu@2336: {
mathieu@2336:   return m_variable->GetValue ();
mathieu@2336: }
mathieu@2336: 
mathieu@2336: uint32_t 
mathieu@2439: RandomVariable::GetInteger (void) const
mathieu@2336: {
mathieu@2439:   return m_variable->GetInteger ();
mathieu@2336: }
raj@4218: 
mathieu@2336: RandomVariableBase *
mathieu@2384: RandomVariable::Peek (void) const
mathieu@2336: {
mathieu@2336:   return m_variable;
mathieu@2336: }
mathieu@2336: 
raj@4218: 
mathieu@2445: ATTRIBUTE_VALUE_IMPLEMENT (RandomVariable);
mathieu@2445: ATTRIBUTE_CHECKER_IMPLEMENT (RandomVariable);
mathieu@2336: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // UniformVariableImpl
mathieu@2336: 
mathieu@2336: class UniformVariableImpl : public RandomVariableBase {
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Creates a uniform random number generator in the
mathieu@2336:    * range [0.0 .. 1.0).
mathieu@2336:    */
mathieu@2336:   UniformVariableImpl();
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Creates a uniform random number generator with the specified range
mathieu@2336:    * \param s Low end of the range
mathieu@2336:    * \param l High end of the range
mathieu@2336:    */
mathieu@2336:   UniformVariableImpl(double s, double l);
mathieu@2336: 
mathieu@2336:   UniformVariableImpl(const UniformVariableImpl& c);
mathieu@2384: 
mathieu@2384:   double GetMin (void) const;
mathieu@2384:   double GetMax (void) const;
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A value between low and high values specified by the constructor
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
raj@4218: 
raj@4218:   /**
raj@4218:    * \return A value between low and high values specified by parameters
raj@4218:    */
raj@4218:   virtual double GetValue(double s, double l);
raj@4218:   
mathieu@2336:   virtual RandomVariableBase*  Copy(void) const;
mathieu@2336: 
mathieu@2336: private:
mathieu@2336:   double m_min;
mathieu@2336:   double m_max;
mathieu@2336: };
mathieu@2336: 
mathieu@2336: UniformVariableImpl::UniformVariableImpl() 
tomh@346:   : m_min(0), m_max(1.0) { }
tomh@346:   
mathieu@2336: UniformVariableImpl::UniformVariableImpl(double s, double l) 
tomh@346:   : m_min(s), m_max(l) { }
tomh@346: 
mathieu@2336: UniformVariableImpl::UniformVariableImpl(const UniformVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max) { }
tomh@346: 
mathieu@2384: double 
mathieu@2384: UniformVariableImpl::GetMin (void) const
mathieu@2384: {
mathieu@2384:   return m_min;
mathieu@2384: }
mathieu@2384: double 
mathieu@2384: UniformVariableImpl::GetMax (void) const
mathieu@2384: {
mathieu@2384:   return m_max;
mathieu@2384: }
mathieu@2384: 
mathieu@2384: 
mathieu@2336: double UniformVariableImpl::GetValue()
tomh@346: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
tomh@346:   return m_min + m_generator->RandU01() * (m_max - m_min);
tomh@346: }
tomh@346: 
raj@4218: double UniformVariableImpl::GetValue(double s, double l) 
raj@4218: {
raj@4224:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
raj@4224:     return s + m_generator->RandU01() * (l-s); 
raj@4218: }
raj@4218: 
mathieu@2336: RandomVariableBase* UniformVariableImpl::Copy() const
mathieu@2336: {
mathieu@2336:   return new UniformVariableImpl(*this);
mathieu@2336: }
mathieu@2336: 
mathieu@2336: UniformVariable::UniformVariable()
mathieu@2336:   : RandomVariable (UniformVariableImpl ())
mathieu@2336: {}
mathieu@2336: UniformVariable::UniformVariable(double s, double l)
mathieu@2336:   : RandomVariable (UniformVariableImpl (s, l))
mathieu@2336: {}
raj@4218: 
mathieu@4242: double UniformVariable::GetValue(void) const
mathieu@2336: {
mathieu@4242:   return this->RandomVariable::GetValue ();
raj@4218: }
raj@4218: 
raj@4218: double UniformVariable::GetValue(double s, double l)
raj@4218: {
raj@4218:   return ((UniformVariableImpl*)Peek())->GetValue(s,l);
tomh@346: }
raj@442: 
timo@4256: uint32_t UniformVariable::GetInteger (uint32_t s, uint32_t l)
timo@4256: {
timo@4256:   NS_ASSERT(s <= l);
timo@4256:   return static_cast<uint32_t>( GetValue(s, l+1) );
timo@4256: }
mathieu@2336: 
mathieu@2336: //-----------------------------------------------------------------------------
mathieu@2336: //-----------------------------------------------------------------------------
mathieu@2336: // ConstantVariableImpl methods
mathieu@2336: 
mathieu@2336: class ConstantVariableImpl : public RandomVariableBase { 
mathieu@2336: 
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Construct a ConstantVariableImpl RNG that returns zero every sample
mathieu@2336:    */
mathieu@2336:   ConstantVariableImpl();
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * Construct a ConstantVariableImpl RNG that returns the specified value
mathieu@2336:    * every sample.
mathieu@2336:    * \param c Unchanging value for this RNG.
mathieu@2336:    */
mathieu@2336:   ConstantVariableImpl(double c);
mathieu@2336: 
mathieu@2336: 
mathieu@2336:   ConstantVariableImpl(const ConstantVariableImpl& c) ;
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \brief Specify a new constant RNG for this generator.
mathieu@2336:    * \param c New constant value for this RNG.
mathieu@2336:    */
mathieu@2336:   void    NewConstant(double c);
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \return The constant value specified
mathieu@2336:    */
mathieu@2336:   virtual double  GetValue();
mathieu@2439:   virtual uint32_t GetInteger();
mathieu@2336:   virtual RandomVariableBase*   Copy(void) const;
mathieu@2336: private:
mathieu@2336:   double m_const;
mathieu@2336: };
mathieu@2336: 
mathieu@2336: ConstantVariableImpl::ConstantVariableImpl() 
mathieu@2336:   : m_const(0) { }
mathieu@2336: 
mathieu@2336: ConstantVariableImpl::ConstantVariableImpl(double c) 
mathieu@2336:   : m_const(c) { };
mathieu@2336:   
mathieu@2336: ConstantVariableImpl::ConstantVariableImpl(const ConstantVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_const(c.m_const) { }
mathieu@2336: 
mathieu@2336: void ConstantVariableImpl::NewConstant(double c) 
mathieu@2336:   { m_const = c;}
mathieu@2336:   
mathieu@2336: double ConstantVariableImpl::GetValue()
raj@442: {
mathieu@2336:   return m_const;
mathieu@2336: }
mathieu@2336: 
mathieu@2439: uint32_t ConstantVariableImpl::GetInteger()
mathieu@2336: {
mathieu@2336:   return (uint32_t)m_const;
mathieu@2336: }
mathieu@2336: 
mathieu@2336: RandomVariableBase* ConstantVariableImpl::Copy() const
mathieu@2336: {
mathieu@2336:   return new ConstantVariableImpl(*this);
mathieu@2336: }
mathieu@2336: 
mathieu@2336: ConstantVariable::ConstantVariable()
mathieu@2336:   : RandomVariable (ConstantVariableImpl ())
mathieu@2336: {}
mathieu@2336: ConstantVariable::ConstantVariable(double c)
mathieu@2336:   : RandomVariable (ConstantVariableImpl (c))
mathieu@2336: {}
mathieu@2336: void 
mathieu@2336: ConstantVariable::SetConstant(double c)
mathieu@2336: {
mathieu@2336:   *this = ConstantVariable (c);
raj@442: }
raj@442: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // SequentialVariableImpl methods
tomh@346: 
mathieu@2336: 
mathieu@2336: class SequentialVariableImpl : public RandomVariableBase {
tomh@346: 
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * \brief Constructor for the SequentialVariableImpl RNG.
mathieu@2336:    *
mathieu@2336:    * The four parameters define the sequence.  For example
mathieu@2336:    * SequentialVariableImpl(0,5,1,2) creates a RNG that has the sequence
mathieu@2336:    * 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0 ...
mathieu@2336:    * \param f First value of the sequence.
mathieu@2336:    * \param l One more than the last value of the sequence.
mathieu@2336:    * \param i Increment between sequence values
mathieu@2336:    * \param c Number of times each member of the sequence is repeated
mathieu@2336:    */
mathieu@2336:   SequentialVariableImpl(double f, double l, double i = 1, uint32_t c = 1);
tomh@346: 
mathieu@2336:   /**
mathieu@2336:    * \brief Constructor for the SequentialVariableImpl RNG.
mathieu@2336:    *
mathieu@2336:    * Differs from the first only in that the increment parameter is a
mathieu@2336:    * random variable
mathieu@2336:    * \param f First value of the sequence.
mathieu@2336:    * \param l One more than the last value of the sequence.
mathieu@2336:    * \param i Reference to a RandomVariableBase for the sequence increment
mathieu@2336:    * \param c Number of times each member of the sequence is repeated
mathieu@2336:    */
mathieu@2336:   SequentialVariableImpl(double f, double l, const RandomVariable& i, uint32_t c = 1);
tomh@346: 
mathieu@2336:   SequentialVariableImpl(const SequentialVariableImpl& c);
mathieu@2336:   
mathieu@2336:   ~SequentialVariableImpl();
mathieu@2336:   /**
mathieu@2336:    * \return The next value in the Sequence
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase*  Copy(void) const;
mathieu@2336: private:
mathieu@2336:   double m_min;
mathieu@2336:   double m_max;
mathieu@2336:   RandomVariable  m_increment;
mathieu@2336:   uint32_t  m_consecutive;
mathieu@2336:   double m_current;
mathieu@2336:   uint32_t  m_currentConsecutive;
mathieu@2336: };
tomh@346: 
mathieu@2336: SequentialVariableImpl::SequentialVariableImpl(double f, double l, double i, uint32_t c)
mathieu@2336:   : m_min(f), m_max(l), m_increment(ConstantVariable(i)), m_consecutive(c),
tomh@346:     m_current(f), m_currentConsecutive(0)
mathieu@2336: {}
tomh@346: 
mathieu@2336: SequentialVariableImpl::SequentialVariableImpl(double f, double l, const RandomVariable& i, uint32_t c)
mathieu@2336:   : m_min(f), m_max(l), m_increment(i), m_consecutive(c),
mathieu@2336:     m_current(f), m_currentConsecutive(0)
mathieu@2336: {}
mathieu@2336: 
mathieu@2336: SequentialVariableImpl::SequentialVariableImpl(const SequentialVariableImpl& c)
mathieu@2336:   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max),
mathieu@2336:     m_increment(c.m_increment), m_consecutive(c.m_consecutive),
tomh@346:     m_current(c.m_current), m_currentConsecutive(c.m_currentConsecutive)
mathieu@2336: {}
tomh@346: 
mathieu@2336: SequentialVariableImpl::~SequentialVariableImpl()
mathieu@2336: {}
raj@360: 
mathieu@2336: double SequentialVariableImpl::GetValue()
tomh@346: { // Return a sequential series of values
tomh@346:   double r = m_current;
tomh@346:   if (++m_currentConsecutive == m_consecutive)
tomh@346:     { // Time to advance to next
tomh@346:       m_currentConsecutive = 0;
mathieu@2336:       m_current += m_increment.GetValue();
tomh@346:       if (m_current >= m_max)
tomh@346:         m_current = m_min + (m_current - m_max);
tomh@346:     }
tomh@346:   return r;
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* SequentialVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new SequentialVariableImpl(*this);
tomh@346: }
mathieu@2336: 
mathieu@2336: SequentialVariable::SequentialVariable(double f, double l, double i, uint32_t c)
mathieu@2336:   : RandomVariable (SequentialVariableImpl (f, l, i, c))
mathieu@2336: {}
mathieu@2336: SequentialVariable::SequentialVariable(double f, double l, const RandomVariable& i, uint32_t c)
mathieu@2336:   : RandomVariable (SequentialVariableImpl (f, l, i, c))
mathieu@2336: {}
mathieu@2336: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // ExponentialVariableImpl methods
mathieu@2336: 
mathieu@2336: class ExponentialVariableImpl : public RandomVariableBase { 
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Constructs an exponential random variable  with a mean
mathieu@2336:    * value of 1.0.
mathieu@2336:    */
mathieu@2336:   ExponentialVariableImpl();
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \brief Constructs an exponential random variable with a specified mean
mathieu@2336:    * \param m Mean value for the random variable
mathieu@2336:    */
mathieu@2336:   explicit ExponentialVariableImpl(double m);
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \brief Constructs an exponential random variable with spefified
mathieu@2336:    * \brief mean and upper limit.
mathieu@2336:    *
mathieu@2336:    * Since exponential distributions can theoretically return unbounded values,
mathieu@2336:    * it is sometimes useful to specify a fixed upper limit.  Note however when
mathieu@2336:    * the upper limit is specified, the true mean of the distribution is 
mathieu@2336:    * slightly smaller than the mean value specified.
mathieu@2336:    * \param m Mean value of the random variable
mathieu@2336:    * \param b Upper bound on returned values
mathieu@2336:    */
mathieu@2336:   ExponentialVariableImpl(double m, double b);
mathieu@2336: 
mathieu@2336:   ExponentialVariableImpl(const ExponentialVariableImpl& c);
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A random value from this exponential distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
raj@4218: 
mathieu@2336: private:
mathieu@2336:   double m_mean;  // Mean value of RV
mathieu@2336:   double m_bound; // Upper bound on value (if non-zero)
mathieu@2336: };
mathieu@2336: 
mathieu@2336: ExponentialVariableImpl::ExponentialVariableImpl() 
tomh@346:   : m_mean(1.0), m_bound(0) { }
tomh@346:   
mathieu@2336: ExponentialVariableImpl::ExponentialVariableImpl(double m) 
tomh@346:   : m_mean(m), m_bound(0) { }
tomh@346:   
mathieu@2336: ExponentialVariableImpl::ExponentialVariableImpl(double m, double b) 
tomh@346:   : m_mean(m), m_bound(b) { }
tomh@346:   
mathieu@2336: ExponentialVariableImpl::ExponentialVariableImpl(const ExponentialVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_mean(c.m_mean), m_bound(c.m_bound) { }
tomh@346: 
mathieu@2336: double ExponentialVariableImpl::GetValue()
tomh@346: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
raj@4227:   while(1)
raj@4224:     {
raj@4227:       double r = -m_mean*log(m_generator->RandU01());
raj@4227:       if (m_bound == 0 || r <= m_bound) return r;
raj@4227:       //otherwise, try again
raj@4224:     }
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* ExponentialVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new ExponentialVariableImpl(*this);
tomh@346: }
mathieu@2336: 
mathieu@2336: ExponentialVariable::ExponentialVariable()
mathieu@2336:   : RandomVariable (ExponentialVariableImpl ())
mathieu@2336: {}
mathieu@2336: ExponentialVariable::ExponentialVariable(double m)
mathieu@2336:   : RandomVariable (ExponentialVariableImpl (m))
mathieu@2336: {}
mathieu@2336: ExponentialVariable::ExponentialVariable(double m, double b)
mathieu@2336:   : RandomVariable (ExponentialVariableImpl (m, b))
mathieu@2336: {}
mathieu@2336: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // ParetoVariableImpl methods
mathieu@2336: class ParetoVariableImpl : public RandomVariableBase {
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Constructs a pareto random variable with a mean of 1 and a shape
mathieu@2336:    * parameter of 1.5
mathieu@2336:    */
mathieu@2336:   ParetoVariableImpl();
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Constructs a pareto random variable with specified mean and shape
mathieu@2336:    * parameter of 1.5
mathieu@2336:    * \param m Mean value of the distribution
mathieu@2336:    */
mathieu@2336:   explicit ParetoVariableImpl(double m);
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Constructs a pareto random variable with the specified mean value and
mathieu@2336:    * shape parameter.
mathieu@2336:    * \param m Mean value of the distribution
mathieu@2336:    * \param s Shape parameter for the distribution
mathieu@2336:    */
mathieu@2336:   ParetoVariableImpl(double m, double s);
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \brief Constructs a pareto random variable with the specified mean
mathieu@2336:    * \brief value, shape (alpha), and upper bound.
mathieu@2336:    *
mathieu@2336:    * Since pareto distributions can theoretically return unbounded values,
mathieu@2336:    * it is sometimes useful to specify a fixed upper limit.  Note however
mathieu@2336:    * when the upper limit is specified, the true mean of the distribution
mathieu@2336:    * is slightly smaller than the mean value specified.
mathieu@2336:    * \param m Mean value
mathieu@2336:    * \param s Shape parameter
mathieu@2336:    * \param b Upper limit on returned values
mathieu@2336:    */
mathieu@2336:   ParetoVariableImpl(double m, double s, double b);
mathieu@2336: 
mathieu@2336:   ParetoVariableImpl(const ParetoVariableImpl& c);
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A random value from this Pareto distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy() const;
raj@4218: 
mathieu@2336: private:
mathieu@2336:   double m_mean;  // Mean value of RV
mathieu@2336:   double m_shape; // Shape parameter
mathieu@2336:   double m_bound; // Upper bound on value (if non-zero)
mathieu@2336: };
mathieu@2336: 
mathieu@2336: ParetoVariableImpl::ParetoVariableImpl() 
tomh@346:   : m_mean(1.0), m_shape(1.5), m_bound(0) { }
tomh@346: 
mathieu@2336: ParetoVariableImpl::ParetoVariableImpl(double m) 
tomh@346:   : m_mean(m), m_shape(1.5), m_bound(0) { }
tomh@346: 
mathieu@2336: ParetoVariableImpl::ParetoVariableImpl(double m, double s) 
tomh@346:     : m_mean(m), m_shape(s), m_bound(0) { }
tomh@346: 
mathieu@2336: ParetoVariableImpl::ParetoVariableImpl(double m, double s, double b) 
tomh@346:   : m_mean(m), m_shape(s), m_bound(b) { }
tomh@346: 
mathieu@2336: ParetoVariableImpl::ParetoVariableImpl(const ParetoVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_mean(c.m_mean), m_shape(c.m_shape), 
raj@360:     m_bound(c.m_bound) { }
tomh@346: 
mathieu@2336: double ParetoVariableImpl::GetValue()
tomh@346: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
tomh@346:   double scale = m_mean * ( m_shape - 1.0) / m_shape;
raj@4227:   while(1)
raj@4224:     {
raj@4227:       double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
raj@4227:       if (m_bound == 0 || r <= m_bound) return r;
raj@4227:       //otherwise, try again
raj@4224:     }
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* ParetoVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new ParetoVariableImpl(*this);
tomh@346: }
raj@442: 
mathieu@2336: ParetoVariable::ParetoVariable ()
mathieu@2336:   : RandomVariable (ParetoVariableImpl ())
mathieu@2336: {}
mathieu@2336: ParetoVariable::ParetoVariable(double m)
mathieu@2336:   : RandomVariable (ParetoVariableImpl (m))
mathieu@2336: {}
mathieu@2336: ParetoVariable::ParetoVariable(double m, double s)
mathieu@2336:   : RandomVariable (ParetoVariableImpl (m, s))
mathieu@2336: {}
mathieu@2336: ParetoVariable::ParetoVariable(double m, double s, double b)
mathieu@2336:   : RandomVariable (ParetoVariableImpl (m, s, b))
mathieu@2336: {}
mathieu@2336: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // WeibullVariableImpl methods
mathieu@2336: 
mathieu@2336: class WeibullVariableImpl : public RandomVariableBase {
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Constructs a weibull random variable  with a mean
mathieu@2336:    * value of 1.0 and a shape (alpha) parameter of 1
mathieu@2336:    */
mathieu@2336:   WeibullVariableImpl();
mathieu@2336: 
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Constructs a weibull random variable with the specified mean
mathieu@2336:    * value and a shape (alpha) parameter of 1.5.
mathieu@2336:    * \param m mean value of the distribution
mathieu@2336:    */
mathieu@2336:    WeibullVariableImpl(double m) ;
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Constructs a weibull random variable with the specified mean
mathieu@2336:    * value and a shape (alpha).
mathieu@2336:    * \param m Mean value for the distribution.
mathieu@2336:    * \param s Shape (alpha) parameter for the distribution.
mathieu@2336:    */
mathieu@2336:   WeibullVariableImpl(double m, double s);
mathieu@2336: 
mathieu@2336:    /**
mathieu@2336:    * \brief Constructs a weibull random variable with the specified mean
mathieu@2336:    * \brief value, shape (alpha), and upper bound.
mathieu@2336:    * Since WeibullVariableImpl distributions can theoretically return unbounded values,
mathieu@2336:    * it is sometimes usefull to specify a fixed upper limit.  Note however
mathieu@2336:    * that when the upper limit is specified, the true mean of the distribution
mathieu@2336:    * is slightly smaller than the mean value specified.
mathieu@2336:    * \param m Mean value for the distribution.
mathieu@2336:    * \param s Shape (alpha) parameter for the distribution.
mathieu@2336:    * \param b Upper limit on returned values
mathieu@2336:    */
mathieu@2336:   WeibullVariableImpl(double m, double s, double b);
mathieu@2336: 
mathieu@2336:   WeibullVariableImpl(const WeibullVariableImpl& c);
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A random value from this Weibull distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
raj@4218: 
mathieu@2336: private:
mathieu@2336:   double m_mean;  // Mean value of RV
mathieu@2336:   double m_alpha; // Shape parameter
mathieu@2336:   double m_bound; // Upper bound on value (if non-zero)
mathieu@2336: };
mathieu@2336: 
mathieu@2336: WeibullVariableImpl::WeibullVariableImpl() : m_mean(1.0), m_alpha(1), m_bound(0) { }
mathieu@2336: WeibullVariableImpl::WeibullVariableImpl(double m) 
tomh@346:   : m_mean(m), m_alpha(1), m_bound(0) { }
mathieu@2336: WeibullVariableImpl::WeibullVariableImpl(double m, double s) 
tomh@346:   : m_mean(m), m_alpha(s), m_bound(0) { }
mathieu@2336: WeibullVariableImpl::WeibullVariableImpl(double m, double s, double b) 
tomh@346:   : m_mean(m), m_alpha(s), m_bound(b) { };
mathieu@2336: WeibullVariableImpl::WeibullVariableImpl(const WeibullVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_mean(c.m_mean), m_alpha(c.m_alpha),
raj@360:     m_bound(c.m_bound) { }
tomh@346: 
mathieu@2336: double WeibullVariableImpl::GetValue()
tomh@346: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
tomh@346:   double exponent = 1.0 / m_alpha;
raj@4227:   while(1)
raj@4224:     {
raj@4227:       double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
raj@4227:       if (m_bound == 0 || r <= m_bound) return r;
raj@4227:       //otherwise, try again
raj@4224:     }
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* WeibullVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new WeibullVariableImpl(*this);
tomh@346: }
raj@442: 
mathieu@2336: WeibullVariable::WeibullVariable()
mathieu@2336:   : RandomVariable (WeibullVariableImpl ())
mathieu@2336: {}
mathieu@2336: WeibullVariable::WeibullVariable(double m)
mathieu@2336:   : RandomVariable (WeibullVariableImpl (m))
mathieu@2336: {}
mathieu@2336: WeibullVariable::WeibullVariable(double m, double s)
mathieu@2336:   : RandomVariable (WeibullVariableImpl (m, s))
mathieu@2336: {}
mathieu@2336: WeibullVariable::WeibullVariable(double m, double s, double b)
mathieu@2336:   : RandomVariable (WeibullVariableImpl (m, s, b))
mathieu@2336: {}
raj@4218: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // NormalVariableImpl methods
mathieu@2336: 
mathieu@2336: class NormalVariableImpl : public RandomVariableBase { // Normally Distributed random var
mathieu@2336: 
mathieu@2336: public:
mathieu@2336:    static const double INFINITE_VALUE;
mathieu@2336:   /**
mathieu@2336:    * Constructs an normal random variable  with a mean
mathieu@2336:    * value of 0 and variance of 1.
mathieu@2336:    */ 
mathieu@2336:   NormalVariableImpl();
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * \brief Construct a normal random variable with specified mean and variance
mathieu@2336:    * \param m Mean value
mathieu@2336:    * \param v Variance
raj@4227:    * \param b Bound.  The NormalVariableImpl is bounded within +-bound of the mean.
mathieu@2336:    */ 
mathieu@2336:   NormalVariableImpl(double m, double v, double b = INFINITE_VALUE);
raj@442: 
mathieu@2336:   NormalVariableImpl(const NormalVariableImpl& c);
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A value from this normal distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
raj@4218: 
gjc@4249:   double GetMean (void) const;
gjc@4249:   double GetVariance (void) const;
gjc@4249:   double GetBound (void) const;
gjc@4249: 
mathieu@2336: private:
mathieu@2336:   double m_mean;      // Mean value of RV
mathieu@2336:   double m_variance;  // Mean value of RV
raj@4227:   double m_bound;     // Bound on value's difference from the mean (absolute value)
mathieu@2336:   bool   m_nextValid; // True if next valid
mathieu@2336:   double m_next;      // The algorithm produces two values at a time
mathieu@2336:   static bool   m_static_nextValid;
mathieu@2336:   static double m_static_next;
mathieu@2336: };
mathieu@2336: 
mathieu@2336: bool         NormalVariableImpl::m_static_nextValid = false;
mathieu@2336: double       NormalVariableImpl::m_static_next;
mathieu@2336: const double NormalVariableImpl::INFINITE_VALUE = 1e307;
mathieu@2336: 
mathieu@2336: NormalVariableImpl::NormalVariableImpl() 
tomh@346:   : m_mean(0.0), m_variance(1.0), m_bound(INFINITE_VALUE), m_nextValid(false){}
tomh@346: 
mathieu@2336: NormalVariableImpl::NormalVariableImpl(double m, double v, double b/*=INFINITE_VALUE*/)
tomh@346:   : m_mean(m), m_variance(v), m_bound(b), m_nextValid(false) { }
tomh@346: 
mathieu@2336: NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c)
mathieu@2336:   : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance),
raj@4218:     m_bound(c.m_bound) { }
tomh@346: 
mathieu@2336: double NormalVariableImpl::GetValue()
tomh@346: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
tomh@346:   if (m_nextValid)
tomh@346:     { // use previously generated
tomh@346:       m_nextValid = false;
tomh@346:       return m_next;
tomh@346:     }
tomh@346:   while(1)
tomh@346:     { // See Simulation Modeling and Analysis p. 466 (Averill Law)
craigdo@3292:       // for algorithm; basically a Box-Muller transform:
craigdo@3292:       // http://en.wikipedia.org/wiki/Box-Muller_transform
tomh@346:       double u1 = m_generator->RandU01();
tomh@346:       double u2 = m_generator->RandU01();;
tomh@346:       double v1 = 2 * u1 - 1;
tomh@346:       double v2 = 2 * u2 - 1;
tomh@346:       double w = v1 * v1 + v2 * v2;
tomh@346:       if (w <= 1.0)
tomh@346:         { // Got good pair
tomh@346:           double y = sqrt((-2 * log(w))/w);
tomh@346:           m_next = m_mean + v2 * y * sqrt(m_variance);
raj@4227:           //if next is in bounds, it is valid
raj@4227:           m_nextValid = fabs(m_next-m_mean) <= m_bound;
tomh@346:           double x1 = m_mean + v1 * y * sqrt(m_variance);
raj@4227:           //if x1 is in bounds, return it
raj@4227:           if (fabs(x1-m_mean) <= m_bound)
raj@4224:             {
raj@4227:               return x1;
raj@4224:             }
raj@4227:           //otherwise try and return m_next if it is valid
raj@4227:           else if (m_nextValid)
raj@4227: 	    {
raj@4227: 	      m_nextValid = false;
raj@4227: 	      return m_next;
raj@4227: 	    }
raj@4227:           //otherwise, just run this loop again
tomh@346:         }
tomh@346:     }
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* NormalVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new NormalVariableImpl(*this);
tomh@346: }
tomh@346: 
gjc@4249: double
gjc@4249: NormalVariableImpl::GetMean (void) const
gjc@4249: {
gjc@4249:   return m_mean;
gjc@4249: }
gjc@4249: 
gjc@4249: double
gjc@4249: NormalVariableImpl::GetVariance (void) const
gjc@4249: {
gjc@4249:   return m_variance;
gjc@4249: }
gjc@4249: 
gjc@4249: double
gjc@4249: NormalVariableImpl::GetBound (void) const
gjc@4249: {
gjc@4249:   return m_bound;
gjc@4249: }
gjc@4249: 
mathieu@2336: NormalVariable::NormalVariable()
mathieu@2336:   : RandomVariable (NormalVariableImpl ())
mathieu@2336: {}
raj@4223: NormalVariable::NormalVariable(double m, double v)
raj@4223:   : RandomVariable (NormalVariableImpl (m, v))
raj@4223: {}
mathieu@2336: NormalVariable::NormalVariable(double m, double v, double b)
mathieu@2336:   : RandomVariable (NormalVariableImpl (m, v, b))
mathieu@2336: {}
mathieu@2336: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: class EmpiricalVariableImpl : public RandomVariableBase {
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Constructor for the EmpiricalVariableImpl random variables.
mathieu@2336:    */
mathieu@2336:   explicit EmpiricalVariableImpl();
mathieu@2336: 
mathieu@2336:   virtual ~EmpiricalVariableImpl();
mathieu@2336:   EmpiricalVariableImpl(const EmpiricalVariableImpl& c);
mathieu@2336:   /**
mathieu@2336:    * \return A value from this empirical distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
mathieu@2336:   /**
mathieu@2336:    * \brief Specifies a point in the empirical distribution
mathieu@2336:    * \param v The function value for this point
mathieu@2336:    * \param c Probability that the function is less than or equal to v
mathieu@2336:    */
mathieu@2336:   virtual void CDF(double v, double c);  // Value, prob <= Value
mathieu@2336: 
mathieu@2336: private:
mathieu@2336:   class ValueCDF {
mathieu@2336:   public:
mathieu@2336:     ValueCDF();
mathieu@2336:     ValueCDF(double v, double c);
mathieu@2336:     ValueCDF(const ValueCDF& c);
mathieu@2336:     double value;
mathieu@2336:     double    cdf;
mathieu@2336:   };
mathieu@2336:   virtual void Validate();  // Insure non-decreasing emiprical values
mathieu@2336:   virtual double Interpolate(double, double, double, double, double);
mathieu@2336:   bool validated; // True if non-decreasing validated
mathieu@2336:   std::vector<ValueCDF> emp;       // Empicical CDF
mathieu@2336: };
mathieu@2336: 
mathieu@2336: 
tomh@346: // ValueCDF methods
mathieu@2336: EmpiricalVariableImpl::ValueCDF::ValueCDF() 
tomh@346:   : value(0.0), cdf(0.0){ }
mathieu@2336: EmpiricalVariableImpl::ValueCDF::ValueCDF(double v, double c) 
tomh@346:   : value(v), cdf(c) { }
mathieu@2336: EmpiricalVariableImpl::ValueCDF::ValueCDF(const ValueCDF& c) 
tomh@346:   : value(c.value), cdf(c.cdf) { }
tomh@346: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2336: // EmpiricalVariableImpl methods
mathieu@2336: EmpiricalVariableImpl::EmpiricalVariableImpl() 
tomh@346:   : validated(false) { }
tomh@346: 
mathieu@2336: EmpiricalVariableImpl::EmpiricalVariableImpl(const EmpiricalVariableImpl& c)
mathieu@2336:   : RandomVariableBase(c), validated(c.validated), emp(c.emp) { }
tomh@346: 
mathieu@2336: EmpiricalVariableImpl::~EmpiricalVariableImpl() { }
tomh@346: 
mathieu@2336: double EmpiricalVariableImpl::GetValue()
tomh@346: { // Return a value from the empirical distribution
tomh@346:   // This code based (loosely) on code by Bruce Mah (Thanks Bruce!)
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
raj@4224:   if (emp.size() == 0) 
raj@4224:     {
raj@4224:       return 0.0; // HuH? No empirical data
raj@4224:     }
raj@4224:   if (!validated) 
raj@4224:     {
raj@4224:       Validate();      // Insure in non-decreasing
raj@4224:     }
tomh@346:   double r = m_generator->RandU01();
raj@4224:   if (r <= emp.front().cdf)
raj@4224:     {
raj@4224:       return emp.front().value; // Less than first
raj@4224:     }
raj@4224:   if (r >= emp.back().cdf) 
raj@4224:     { 
raj@4224:       return emp.back().value;  // Greater than last
raj@4224:     }
tomh@346:   // Binary search
tomh@346:   std::vector<ValueCDF>::size_type bottom = 0;
tomh@346:   std::vector<ValueCDF>::size_type top = emp.size() - 1;
tomh@346:   while(1)
tomh@346:     {
tomh@346:       std::vector<ValueCDF>::size_type c = (top + bottom) / 2;
tomh@346:       if (r >= emp[c].cdf && r < emp[c+1].cdf)
tomh@346:         { // Found it
tomh@346:           return Interpolate(emp[c].cdf, emp[c+1].cdf,
tomh@346:                              emp[c].value, emp[c+1].value,
tomh@346:                              r);
tomh@346:         }
tomh@346:       // Not here, adjust bounds
raj@4224:       if (r < emp[c].cdf)
raj@4224:         {
raj@4224:           top    = c - 1;
raj@4224:         }
raj@4224:       else
raj@4224:         {
raj@4224:           bottom = c + 1;
raj@4224:         }
tomh@346:     }
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* EmpiricalVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new EmpiricalVariableImpl(*this);
tomh@346: }
tomh@346: 
mathieu@2336: void EmpiricalVariableImpl::CDF(double v, double c)
tomh@346: { // Add a new empirical datapoint to the empirical cdf
tomh@346:   // NOTE.   These MUST be inserted in non-decreasing order
tomh@346:   emp.push_back(ValueCDF(v, c));
tomh@346: }
tomh@346: 
mathieu@2336: void EmpiricalVariableImpl::Validate()
tomh@346: {
tomh@346:   ValueCDF prior;
tomh@346:   for (std::vector<ValueCDF>::size_type i = 0; i < emp.size(); ++i)
tomh@346:     {
tomh@346:       ValueCDF& current = emp[i];
tomh@346:       if (current.value < prior.value || current.cdf < prior.cdf)
tomh@346:         { // Error
gjc@1020:           cerr << "Empirical Dist error,"
tomh@346:                << " current value " << current.value
tomh@346:                << " prior value "   << prior.value
tomh@346:                << " current cdf "   << current.cdf
tomh@346:                << " prior cdf "     << prior.cdf << endl;
tomh@346:           NS_FATAL_ERROR("Empirical Dist error");
tomh@346:         }
tomh@346:       prior = current;
tomh@346:     }
tomh@346:   validated = true;
tomh@346: }
tomh@346: 
mathieu@2336: double EmpiricalVariableImpl::Interpolate(double c1, double c2,
tomh@346:                                 double v1, double v2, double r)
tomh@346: { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)
tomh@346:   return (v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
tomh@346: }
tomh@346: 
mathieu@2336: EmpiricalVariable::EmpiricalVariable()
mathieu@2336:   : RandomVariable (EmpiricalVariableImpl ())
mathieu@2336: {}
mathieu@2336: EmpiricalVariable::EmpiricalVariable (const RandomVariableBase &variable)
mathieu@2336:   : RandomVariable (variable)
mathieu@2336: {}
mathieu@2336: void 
mathieu@2336: EmpiricalVariable::CDF(double v, double c)
tomh@346: {
mathieu@2336:   EmpiricalVariableImpl *impl = dynamic_cast<EmpiricalVariableImpl *> (Peek ());
mathieu@2336:   NS_ASSERT (impl);
mathieu@2336:   impl->CDF (v, c);
tomh@346: }
tomh@346: 
tomh@346: 
tomh@346: //-----------------------------------------------------------------------------
tomh@346: //-----------------------------------------------------------------------------
mathieu@2965: // IntegerValue EmpiricalVariableImpl methods
mathieu@2336: class IntEmpiricalVariableImpl : public EmpiricalVariableImpl {
mathieu@2336: public:
mathieu@2336: 
mathieu@2336:   IntEmpiricalVariableImpl();
mathieu@2336:   
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
mathieu@2336:   /**
mathieu@2336:    * \return An integer value from this empirical distribution
mathieu@2336:    */
mathieu@2439:   virtual uint32_t GetInteger();
mathieu@2336: private:
mathieu@2336:   virtual double Interpolate(double, double, double, double, double);
mathieu@2336: };
mathieu@2336: 
mathieu@2336: 
mathieu@2336: IntEmpiricalVariableImpl::IntEmpiricalVariableImpl() { }
mathieu@2336: 
mathieu@2439: uint32_t IntEmpiricalVariableImpl::GetInteger()
mathieu@2336: {
mathieu@2336:   return (uint32_t)GetValue();
mathieu@2336: }
mathieu@2336: 
mathieu@2336: RandomVariableBase* IntEmpiricalVariableImpl::Copy() const
mathieu@2336: {
mathieu@2336:   return new IntEmpiricalVariableImpl(*this);
mathieu@2336: }
mathieu@2336: 
mathieu@2336: double IntEmpiricalVariableImpl::Interpolate(double c1, double c2,
mathieu@2336:                                    double v1, double v2, double r)
mathieu@2336: { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)
mathieu@2336:   return ceil(v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
mathieu@2336: }
mathieu@2336: 
mathieu@2336: IntEmpiricalVariable::IntEmpiricalVariable()
mathieu@2336:   : EmpiricalVariable (IntEmpiricalVariableImpl ())
mathieu@2336: {}
mathieu@2336: 
mathieu@2336: //-----------------------------------------------------------------------------
mathieu@2336: //-----------------------------------------------------------------------------
mathieu@2336: // DeterministicVariableImpl
mathieu@2336: class DeterministicVariableImpl : public RandomVariableBase 
mathieu@2336: {
mathieu@2336: 
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * \brief Constructor
mathieu@2336:    *
mathieu@2336:    * Creates a generator that returns successive elements of the d array
mathieu@2336:    * on successive calls to ::Value().  Note that the d pointer is copied
mathieu@2336:    * for use by the generator (shallow-copy), not its contents, so the 
mathieu@2336:    * contents of the array d points to have to remain unchanged for the use 
mathieu@2336:    * of DeterministicVariableImpl to be meaningful.
mathieu@2336:    * \param d Pointer to array of random values to return in sequence
mathieu@2336:    * \param c Number of values in the array
mathieu@2336:    */
mathieu@2336:   explicit DeterministicVariableImpl(double* d, uint32_t c);
mathieu@2336: 
mathieu@2336:   virtual ~DeterministicVariableImpl();
mathieu@2336:   /**
mathieu@2336:    * \return The next value in the deterministic sequence
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
mathieu@2336: private:
mathieu@2336:   uint32_t   count;
mathieu@2336:   uint32_t   next;
mathieu@2336:   double* data;
mathieu@2336: };
mathieu@2336: 
mathieu@2336: DeterministicVariableImpl::DeterministicVariableImpl(double* d, uint32_t c)
tomh@346:     : count(c), next(c), data(d)
tomh@346: { // Nothing else needed
tomh@346: }
tomh@346: 
mathieu@2336: DeterministicVariableImpl::~DeterministicVariableImpl() { }
tomh@346:   
mathieu@2336: double DeterministicVariableImpl::GetValue()
tomh@346: {
raj@4224:   if (next == count) 
raj@4224:     {
raj@4224:       next = 0;
raj@4224:     }
tomh@346:   return data[next++];
tomh@346: }
tomh@346: 
mathieu@2336: RandomVariableBase* DeterministicVariableImpl::Copy() const
tomh@346: {
mathieu@2336:   return new DeterministicVariableImpl(*this);
tomh@346: }
tomh@346: 
mathieu@2336: DeterministicVariable::DeterministicVariable(double* d, uint32_t c)
mathieu@2336:   : RandomVariable (DeterministicVariableImpl (d, c))
mathieu@2336: {}
raj@438: 
raj@438: //-----------------------------------------------------------------------------
raj@438: //-----------------------------------------------------------------------------
mathieu@2336: // LogNormalVariableImpl
mathieu@2336: class LogNormalVariableImpl : public RandomVariableBase { 
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * \param mu mu parameter of the lognormal distribution
mathieu@2336:    * \param sigma sigma parameter of the lognormal distribution
mathieu@2336:    */
mathieu@2336:   LogNormalVariableImpl (double mu, double sigma);
raj@438: 
mathieu@2336:   /**
mathieu@2336:    * \return A random value from this distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue ();
mathieu@2336:   virtual RandomVariableBase* Copy(void) const;
raj@4218: 
mathieu@2336: private:
mathieu@2336:   double m_mu;
mathieu@2336:   double m_sigma;
mathieu@2336: };
mathieu@2336: 
mathieu@2336: 
mathieu@2336: RandomVariableBase* LogNormalVariableImpl::Copy () const
raj@438: {
francesco@4579:   return new LogNormalVariableImpl (*this);
raj@438: }
raj@438: 
mathieu@2336: LogNormalVariableImpl::LogNormalVariableImpl (double mu, double sigma)
raj@442:     :m_mu(mu), m_sigma(sigma) 
raj@438: {
raj@438: }
raj@438: 
raj@438: // The code from this function was adapted from the GNU Scientific
raj@438: // Library 1.8:
raj@438: /* randist/lognormal.c
raj@438:  * 
raj@438:  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
raj@438:  * 
raj@438:  * This program is free software; you can redistribute it and/or modify
raj@438:  * it under the terms of the GNU General Public License as published by
raj@438:  * the Free Software Foundation; either version 2 of the License, or (at
raj@438:  * your option) any later version.
raj@438:  * 
raj@438:  * This program is distributed in the hope that it will be useful, but
raj@438:  * WITHOUT ANY WARRANTY; without even the implied warranty of
raj@438:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
raj@438:  * General Public License for more details.
raj@438:  * 
raj@438:  * You should have received a copy of the GNU General Public License
raj@438:  * along with this program; if not, write to the Free Software
raj@438:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
raj@438:  */
raj@438: /* The lognormal distribution has the form 
raj@438: 
raj@438:    p(x) dx = 1/(x * sqrt(2 pi sigma^2)) exp(-(ln(x) - zeta)^2/2 sigma^2) dx
raj@438: 
raj@438:    for x > 0. Lognormal random numbers are the exponentials of
raj@438:    gaussian random numbers */
raj@438: double
mathieu@2336: LogNormalVariableImpl::GetValue ()
raj@438: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
raj@438:   double u, v, r2, normal, z;
raj@438: 
raj@438:   do
raj@438:     {
raj@438:       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
raj@438: 
raj@438:       u = -1 + 2 * m_generator->RandU01 ();
raj@438:       v = -1 + 2 * m_generator->RandU01 ();
raj@438: 
raj@438:       /* see if it is in the unit circle */
raj@438:       r2 = u * u + v * v;
raj@438:     }
raj@438:   while (r2 > 1.0 || r2 == 0);
raj@438: 
raj@438:   normal = u * sqrt (-2.0 * log (r2) / r2);
raj@438: 
raj@438:   z =  exp (m_sigma * normal + m_mu);
raj@438: 
raj@438:   return z;
raj@438: }
raj@438: 
mathieu@2336: LogNormalVariable::LogNormalVariable (double mu, double sigma)
mathieu@2336:   : RandomVariable (LogNormalVariableImpl (mu, sigma))
mathieu@2336: {}
mathieu@2336: 
raj@942: //-----------------------------------------------------------------------------
raj@942: //-----------------------------------------------------------------------------
timo@4317: // GammaVariableImpl
timo@4317: class GammaVariableImpl : public RandomVariableBase
timo@4317: {
timo@4317: public:
timo@4317:   /**
timo@4317:    * \param alpha alpha parameter of the gamma distribution
timo@4317:    * \param beta beta parameter of the gamma distribution
timo@4317:    */
timo@4317:   GammaVariableImpl (double alpha, double beta);
timo@4317: 
timo@4317:   /**
timo@4317:    * \return A random value from this distribution
timo@4317:    */
timo@4317:   virtual double GetValue ();
timo@4317: 
timo@4317:   /**
timo@4317:    * \return A random value from the gamma distribution with parameters alpha
timo@4317:    * and beta
timo@4317:    */
timo@4317:   double GetValue(double alpha, double beta);
timo@4317: 
timo@4317:   virtual RandomVariableBase* Copy(void) const;
timo@4317: 
timo@4317: private:
timo@4317:   double m_alpha;
timo@4317:   double m_beta;
timo@4317:   NormalVariable m_normal;
timo@4317: };
timo@4317: 
timo@4317: 
timo@4317: RandomVariableBase* GammaVariableImpl::Copy () const
timo@4317: {
timo@4317:   return new GammaVariableImpl (m_alpha, m_beta);
timo@4317: }
timo@4317: 
timo@4317: GammaVariableImpl::GammaVariableImpl (double alpha, double beta)
timo@4317:   : m_alpha(alpha), m_beta(beta) 
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: double
timo@4317: GammaVariableImpl::GetValue ()
timo@4317: {
timo@4317:   return GetValue(m_alpha, m_beta);
timo@4317: }
timo@4317: 
timo@4317: /*
timo@4317:   The code for the following generator functions was adapted from ns-2
timo@4317:   tools/ranvar.cc
timo@4317: 
timo@4317:   Originally the algorithm was devised by Marsaglia in 2000:
timo@4317:   G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
timo@4317:   ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
timo@4317: 
timo@4317:   The Gamma distribution density function has the form 
timo@4317: 
timo@4317: 	                     x^(alpha-1) * exp(-x/beta)
timo@4317: 	p(x; alpha, beta) = ----------------------------
timo@4317:                              beta^alpha * Gamma(alpha)
timo@4317: 
timo@4317:   for x > 0.
timo@4317: */
timo@4317: double
timo@4317: GammaVariableImpl::GetValue (double alpha, double beta)
timo@4317: {
timo@4317:   if(!m_generator)
timo@4317:     {
timo@4317:       m_generator = new RngStream();
timo@4317:     }
timo@4317: 
timo@4317:   if (alpha < 1)
timo@4317:     {
timo@4317:       double u = m_generator->RandU01 ();
timo@4317:       return GetValue(1.0 + alpha, beta) * pow (u, 1.0 / alpha);
timo@4317:     }
timo@4317: 	
timo@4317:   double x, v, u;
timo@4317:   double d = alpha - 1.0 / 3.0;
timo@4317:   double c = (1.0 / 3.0) / sqrt (d);
timo@4317: 
timo@4317:   while (1)
timo@4317:     {
timo@4317:       do
timo@4317:         {
timo@4317:           x = m_normal.GetValue ();
timo@4317:           v = 1.0 + c * x;
timo@4317:         } while (v <= 0);
timo@4317: 
timo@4317:       v = v * v * v;
timo@4317:       u = m_generator->RandU01 ();
timo@4317:       if (u < 1 - 0.0331 * x * x * x * x)
timo@4317:         {
timo@4317:           break;
timo@4317:         }
timo@4317:       if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
timo@4317:         {
timo@4317:           break;
timo@4317:         }
timo@4317:     }
timo@4317: 
timo@4317:   return beta * d * v;
timo@4317: }
timo@4317: 
timo@4317: GammaVariable::GammaVariable ()
timo@4317:   : RandomVariable (GammaVariableImpl (1.0, 1.0))
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: GammaVariable::GammaVariable (double alpha, double beta)
timo@4317:   : RandomVariable (GammaVariableImpl (alpha, beta))
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: double GammaVariable::GetValue(void) const
timo@4317: {
timo@4317:   return this->RandomVariable::GetValue ();
timo@4317: }
timo@4317: 
timo@4317: double GammaVariable::GetValue(double alpha, double beta) const
timo@4317: {
timo@4317:   return ((GammaVariableImpl*)Peek())->GetValue(alpha, beta);
timo@4317: }
timo@4317: 
timo@4317: //-----------------------------------------------------------------------------
timo@4317: //-----------------------------------------------------------------------------
timo@4317: // ErlangVariableImpl
timo@4317: 
timo@4317: class ErlangVariableImpl : public RandomVariableBase
timo@4317: {
timo@4317: public:
timo@4317:   /**
timo@4317:    * \param k k parameter of the Erlang distribution
timo@4317:    * \param lambda lambda parameter of the Erlang distribution
timo@4317:    */
timo@4317:   ErlangVariableImpl (unsigned int k, double lambda);
timo@4317: 
timo@4317:   /**
timo@4317:    * \return A random value from this distribution
timo@4317:    */
timo@4317:   virtual double GetValue ();
timo@4317: 
timo@4317:   /**
timo@4317:    * \return A random value from the Erlang distribution with parameters k and
timo@4317:    * lambda.
timo@4317:    */
timo@4317:   double GetValue(unsigned int k, double lambda);
timo@4317: 
timo@4317:   virtual RandomVariableBase* Copy(void) const;
timo@4317: 
timo@4317: private:
timo@4317:   unsigned int m_k;
timo@4317:   double m_lambda;
timo@4317: };
timo@4317: 
timo@4317: 
timo@4317: RandomVariableBase* ErlangVariableImpl::Copy () const
timo@4317: {
timo@4317:   return new ErlangVariableImpl (m_k, m_lambda);
timo@4317: }
timo@4317: 
timo@4317: ErlangVariableImpl::ErlangVariableImpl (unsigned int k, double lambda)
timo@4317:   : m_k(k), m_lambda(lambda)
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: double
timo@4317: ErlangVariableImpl::GetValue ()
timo@4317: {
timo@4317:   return GetValue(m_k, m_lambda);
timo@4317: }
timo@4317: 
timo@4317: /*
timo@4317:   The code for the following generator functions was adapted from ns-2
timo@4317:   tools/ranvar.cc
timo@4317: 
timo@4317:   The Erlang distribution density function has the form 
timo@4317: 
timo@4317: 	                   x^(k-1) * exp(-x/lambda)
timo@4317: 	p(x; k, lambda) = ---------------------------
timo@4317:                              lambda^k * (k-1)!
timo@4317: 
timo@4317:   for x > 0.
timo@4317: */
timo@4317: double
timo@4317: ErlangVariableImpl::GetValue (unsigned int k, double lambda)
timo@4317: {
timo@4317:   if(!m_generator)
timo@4317:     {
timo@4317:       m_generator = new RngStream();
timo@4317:     }
timo@4317: 
timo@4317:   ExponentialVariable exponential(lambda);
timo@4317: 
timo@4317:   double result = 0;
timo@4317:   for (unsigned int i = 0; i < k; ++i)
timo@4317:     {
timo@4317:       result += exponential.GetValue();
timo@4317:     }
timo@4317: 
timo@4317:   return result;
timo@4317: }
timo@4317: 
timo@4317: ErlangVariable::ErlangVariable ()
timo@4317:   : RandomVariable (ErlangVariableImpl (1, 1.0))
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: ErlangVariable::ErlangVariable (unsigned int k, double lambda)
timo@4317:   : RandomVariable (ErlangVariableImpl (k, lambda))
timo@4317: {
timo@4317: }
timo@4317: 
timo@4317: double ErlangVariable::GetValue(void) const
timo@4317: {
timo@4317:   return this->RandomVariable::GetValue ();
timo@4317: }
timo@4317: 
timo@4317: double ErlangVariable::GetValue(unsigned int k, double lambda) const
timo@4317: {
timo@4317:   return ((ErlangVariableImpl*)Peek())->GetValue(k, lambda);
timo@4317: }
timo@4317: 
timo@4317: //-----------------------------------------------------------------------------
timo@4317: //-----------------------------------------------------------------------------
mathieu@2336: // TriangularVariableImpl methods
mathieu@2336: class TriangularVariableImpl : public RandomVariableBase {
mathieu@2336: public:
mathieu@2336:   /**
mathieu@2336:    * Creates a triangle distribution random number generator in the
mathieu@2336:    * range [0.0 .. 1.0), with mean of 0.5
mathieu@2336:    */
mathieu@2336:   TriangularVariableImpl();
mathieu@2336: 
mathieu@2336:   /**
mathieu@2336:    * Creates a triangle distribution random number generator with the specified
mathieu@2336:    * range
mathieu@2336:    * \param s Low end of the range
mathieu@2336:    * \param l High end of the range
mathieu@2336:    * \param mean mean of the distribution
mathieu@2336:    */
mathieu@2336:   TriangularVariableImpl(double s, double l, double mean);
mathieu@2336: 
mathieu@2336:   TriangularVariableImpl(const TriangularVariableImpl& c);
mathieu@2336:   
mathieu@2336:   /**
mathieu@2336:    * \return A value from this distribution
mathieu@2336:    */
mathieu@2336:   virtual double GetValue();
mathieu@2336:   virtual RandomVariableBase*  Copy(void) const;
raj@4218: 
mathieu@2336: private:
mathieu@2336:   double m_min;
mathieu@2336:   double m_max;
mathieu@2336:   double m_mode;  //easier to work with the mode internally instead of the mean
mathieu@2336:                   //they are related by the simple: mean = (min+max+mode)/3
mathieu@2336: };
mathieu@2336: 
mathieu@2336: TriangularVariableImpl::TriangularVariableImpl() 
raj@942:   : m_min(0), m_max(1), m_mode(0.5) { }
raj@942:   
mathieu@2336: TriangularVariableImpl::TriangularVariableImpl(double s, double l, double mean) 
raj@942:   : m_min(s), m_max(l), m_mode(3.0*mean-s-l) { }
raj@942:   
mathieu@2336: TriangularVariableImpl::TriangularVariableImpl(const TriangularVariableImpl& c) 
mathieu@2336:   : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max), m_mode(c.m_mode) { }
raj@942: 
mathieu@2336: double TriangularVariableImpl::GetValue()
raj@942: {
raj@1806:   if(!m_generator)
raj@4224:     {
raj@4224:       m_generator = new RngStream();
raj@4224:     }
raj@942:   double u = m_generator->RandU01();
raj@942:   if(u <= (m_mode - m_min) / (m_max - m_min) )
raj@4224:     {
raj@4224:       return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) );
raj@4224:     }
raj@942:   else
raj@4224:     {
raj@4224:       return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) );
raj@4224:     }
raj@942: }
raj@942: 
mathieu@2336: RandomVariableBase* TriangularVariableImpl::Copy() const
raj@942: {
mathieu@2336:   return new TriangularVariableImpl(*this);
raj@942: }
raj@942: 
mathieu@2336: TriangularVariable::TriangularVariable()
mathieu@2336:   : RandomVariable (TriangularVariableImpl ())
mathieu@2336: {}
mathieu@2336: TriangularVariable::TriangularVariable(double s, double l, double mean)
mathieu@2336:   : RandomVariable (TriangularVariableImpl (s,l,mean))
mathieu@2336: {}
mathieu@2336: 
francesco@4579: //-----------------------------------------------------------------------------
francesco@4579: //-----------------------------------------------------------------------------
francesco@4579: // ZipfVariableImpl
francesco@4579: class ZipfVariableImpl : public RandomVariableBase { 
francesco@4579: public:
francesco@4579:   /**
francesco@4579:    * \param n the number of possible items
francesco@4579:    * \param alpha the alpha parameter
francesco@4579:    */
francesco@4579:   ZipfVariableImpl (long n, double alpha);
francesco@4579: 
francesco@4579:   /**
francesco@4579:    * \A zipf variable with N=1 and alpha=0
francesco@4579:    */
francesco@4579:   ZipfVariableImpl ();
francesco@4579: 
francesco@4579:   /**
francesco@4579:    * \return A random value from this distribution
francesco@4579:    */
francesco@4579:   virtual double GetValue ();
francesco@4579:   virtual RandomVariableBase* Copy(void) const;
francesco@4579: 
francesco@4579: private:
francesco@4579:   long m_n;
francesco@4579:   double m_alpha;
francesco@4579:   double m_c; //the normalization constant
francesco@4579: };
francesco@4579: 
francesco@4579: 
francesco@4579: RandomVariableBase* ZipfVariableImpl::Copy () const
francesco@4579: {
francesco@4579:   return new ZipfVariableImpl (m_n, m_alpha);
francesco@4579: }
francesco@4579: 
francesco@4579: ZipfVariableImpl::ZipfVariableImpl ()
francesco@4579:     :m_n(1), m_alpha(0), m_c(1)
francesco@4579: {
francesco@4579: }
francesco@4579: 
francesco@4579: 
francesco@4579: ZipfVariableImpl::ZipfVariableImpl (long n, double alpha)
francesco@4579:     :m_n(n), m_alpha(alpha), m_c(0)
francesco@4579: {
francesco@4579:   //calculate the normalization constant c
francesco@4579:   for(int i=1;i<=n;i++)
francesco@4579:     {
francesco@4579:       m_c+=(1.0/pow((double)i,alpha));
francesco@4579:     }
francesco@4579:   m_c=1.0/m_c;
francesco@4579: }
francesco@4579: 
francesco@4579: double
francesco@4579: ZipfVariableImpl::GetValue ()
francesco@4579: {
francesco@4579:   if(!m_generator)
francesco@4579:     {
francesco@4579:       m_generator = new RngStream();
francesco@4579:     }
francesco@4579: 
francesco@4579:   double u = m_generator->RandU01();
francesco@4579:   double sum_prob=0,zipf_value=0;
francesco@4579:   for(int i=1;i<=m_n;i++)
francesco@4579:     {
francesco@4579:       sum_prob+=m_c/pow((double)i,m_alpha);
francesco@4579:       if(sum_prob>u)
francesco@4579:         {
francesco@4579:           zipf_value=i;
francesco@4579:           break;
francesco@4579:         }
francesco@4579:     }
francesco@4579:   return zipf_value;
francesco@4579: }
francesco@4579: 
francesco@4579: ZipfVariable::ZipfVariable ()
francesco@4579:   : RandomVariable (ZipfVariableImpl ())
francesco@4579: {}
francesco@4579: 
francesco@4579: ZipfVariable::ZipfVariable (long n, double alpha)
francesco@4579:   : RandomVariable (ZipfVariableImpl (n, alpha))
francesco@4579: {}
francesco@4579: 
mathieu@2336: 
mathieu@2409: std::ostream &operator << (std::ostream &os, const RandomVariable &var)
mathieu@2384: {
mathieu@2409:   RandomVariableBase *base = var.Peek ();
mathieu@2384:   ConstantVariableImpl *constant = dynamic_cast<ConstantVariableImpl *> (base);
mathieu@2384:   if (constant != 0)
mathieu@2384:     {
mathieu@2409:       os << "Constant:" << constant->GetValue ();
mathieu@2409:       return os;
mathieu@2384:     }
mathieu@2384:   UniformVariableImpl *uniform = dynamic_cast<UniformVariableImpl *> (base);
mathieu@2384:   if (uniform != 0)
mathieu@2384:     {
mathieu@2409:       os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax ();
mathieu@2409:       return os;
mathieu@2384:     }
gjc@4249:   NormalVariableImpl *normal = dynamic_cast<NormalVariableImpl *> (base);
gjc@4249:   if (normal != 0)
gjc@4249:     {
gjc@4249:       os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance ();
gjc@4249:       double bound = normal->GetBound ();
gjc@4249:       if (bound != NormalVariableImpl::INFINITE_VALUE)
gjc@4249:         {
gjc@4249:           os << ":" << bound;
gjc@4249:         }
gjc@4249:       return os;
gjc@4249:     }
mathieu@2384:   // XXX: support other distributions
mathieu@2409:   os.setstate (std::ios_base::badbit);
mathieu@2409:   return os;
mathieu@2384: }
mathieu@2409: std::istream &operator >> (std::istream &is, RandomVariable &var)
mathieu@2384: {
mathieu@2409:   std::string value;
mathieu@2409:   is >> value;
mathieu@2384:   std::string::size_type tmp;
mathieu@2384:   tmp = value.find (":");
mathieu@2384:   if (tmp == std::string::npos)
mathieu@2384:     {
mathieu@2409:       is.setstate (std::ios_base::badbit);
mathieu@2409:       return is;
mathieu@2384:     }
mathieu@2384:   std::string type = value.substr (0, tmp);
gjc@2628:   value = value.substr (tmp + 1, value.npos);
gjc@2628:   if (type == "Constant")
mathieu@2384:     {
gjc@2628:       istringstream iss (value);
gjc@2628:       double constant;
gjc@2628:       iss >> constant;
gjc@2628:       var = ConstantVariable (constant);
mathieu@2384:     }
gjc@2628:   else if (type == "Uniform")
mathieu@2384:     {
gjc@2628:       if (value.size () == 0)
gjc@2628:         {
gjc@2628:           var = UniformVariable ();
gjc@2628:         }
gjc@2628:       else
gjc@2628:         {
gjc@2628:           tmp = value.find (":");
gjc@2628:           if (tmp == value.npos)
gjc@2628:             {
gjc@2628:               NS_FATAL_ERROR ("bad Uniform value: " << value);
gjc@2628:             }
gjc@2628:           istringstream issA (value.substr (0, tmp));
gjc@2628:           istringstream issB (value.substr (tmp + 1, value.npos));
gjc@2628:           double a, b;
gjc@2628:           issA >> a;
gjc@2628:           issB >> b;
gjc@2628:           var = UniformVariable (a, b);
gjc@2628:         }
mathieu@2384:     }
gjc@4249:   else if (type == "Normal")
gjc@4249:     {
gjc@4249:       if (value.size () == 0)
gjc@4249:         {
gjc@4249:           var = NormalVariable ();
gjc@4249:         }
gjc@4249:       else
gjc@4249:         {
gjc@4249:           tmp = value.find (":");
gjc@4249:           if (tmp == value.npos)
gjc@4249:             {
gjc@4249:               NS_FATAL_ERROR ("bad Normal value: " << value);
gjc@4249:             }
gjc@4249:           std::string::size_type tmp2;
gjc@4249:           std::string sub = value.substr (tmp + 1, value.npos);
gjc@4249:           tmp2 = sub.find (":");
gjc@4249:           if (tmp2 == value.npos)
gjc@4249:             {
gjc@4249:               istringstream issA (value.substr (0, tmp));
gjc@4249:               istringstream issB (sub);
gjc@4249:               double a, b;
gjc@4249:               issA >> a;
gjc@4249:               issB >> b;
gjc@4249:               var = NormalVariable (a, b);
gjc@4249:             }
gjc@4249:           else
gjc@4249:             {
gjc@4249:               istringstream issA (value.substr (0, tmp));
gjc@4249:               istringstream issB (sub.substr (0, tmp2));
gjc@4249:               istringstream issC (sub.substr (tmp2 + 1, value.npos));
gjc@4249:               double a, b, c;
gjc@4249:               issA >> a;
gjc@4249:               issB >> b;
gjc@4249:               issC >> c;
gjc@4249:               var = NormalVariable (a, b, c);
gjc@4249:             }
gjc@4249:         }
gjc@4249:     }
mathieu@2409:   else
mathieu@2409:     {
gjc@2628:       NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << type);
mathieu@2409:       // XXX: support other distributions.
mathieu@2409:     }
mathieu@2409:   return is;
mathieu@2384: }
mathieu@2384: 
mathieu@2384: 
raj@942: 
tomh@346: }//namespace ns3
tomh@346: 
mathieu@4243:  
mathieu@4243: 
mathieu@4243: #ifdef RUN_SELF_TESTS
mathieu@4243: #include "test.h"
mathieu@4243: #include <vector>
mathieu@4243: 
mathieu@4243: namespace ns3 {
mathieu@4243: 
mathieu@4243: 
mathieu@4243: class RandomVariableTest : public Test
mathieu@4243: {
mathieu@4243: public:
mathieu@4243:   RandomVariableTest () : Test ("RandomVariable") {}
mathieu@4243:   virtual bool RunTests (void)
mathieu@4243:   {
mathieu@4243:     bool result = true;
mathieu@4243:     const double desired_mean = 1.0;
mathieu@4243:     const double desired_stddev = 1.0;
mathieu@4243:     double tmp = log (1 + (desired_stddev/desired_mean)*(desired_stddev/desired_mean));
mathieu@4243:     double sigma = sqrt (tmp);
mathieu@4243:     double mu = log (desired_mean) - 0.5*tmp;
mathieu@4243: 
mathieu@4243:     // Test a custom lognormal instance
mathieu@4243:     {
mathieu@4243:       LogNormalVariable lognormal (mu, sigma);
mathieu@4243:       vector<double> samples;
mathieu@4243:       const int NSAMPLES = 10000;
mathieu@4243:       double sum = 0;
mathieu@4243:       for (int n = NSAMPLES; n; --n)
mathieu@4243:         {
mathieu@4243:           double value = lognormal.GetValue ();
mathieu@4243:           sum += value;
mathieu@4243:           samples.push_back (value);
mathieu@4243:         }
mathieu@4243:       double obtained_mean = sum / NSAMPLES;
mathieu@4243:       sum = 0;
mathieu@4243:       for (vector<double>::iterator iter = samples.begin (); iter != samples.end (); iter++)
mathieu@4243:         {
mathieu@4243:           double tmp = (*iter - obtained_mean);
mathieu@4243:           sum += tmp*tmp;
mathieu@4243:         }
mathieu@4243:       double obtained_stddev = sqrt (sum / (NSAMPLES - 1));
mathieu@4243: 
mathieu@4243:       if (not (obtained_mean/desired_mean > 0.90 and obtained_mean/desired_mean < 1.10))
mathieu@4243:         {
mathieu@4243:           result = false;
mathieu@4243:           Failure () << "Obtained lognormal mean value " << obtained_mean << ", expected " << desired_mean << std::endl;
mathieu@4243:         }
mathieu@4243: 
mathieu@4243:       if (not (obtained_stddev/desired_stddev > 0.90 and obtained_stddev/desired_stddev < 1.10))
mathieu@4243:         {
mathieu@4243:           result = false;
mathieu@4243:           Failure () << "Obtained lognormal stddev value " << obtained_stddev <<
mathieu@4243:             ", expected " << desired_stddev << std::endl;
mathieu@4243:         }
mathieu@4243:     }
mathieu@4243: 
mathieu@4243:     // Test attribute serialization
mathieu@4243:     {
gjc@4249:       {
gjc@4249:         RandomVariableValue val;
gjc@4249:         val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ());
gjc@4249:         RandomVariable rng = val.Get ();
gjc@4249:         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2");
gjc@4249:       }
gjc@4249:       {
gjc@4249:         RandomVariableValue val;
gjc@4249:         val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ());
gjc@4249:         RandomVariable rng = val.Get ();
gjc@4249:         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2");
gjc@4249:       }
gjc@4249:       {
gjc@4249:         RandomVariableValue val;
gjc@4249:         val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker ());
gjc@4249:         RandomVariable rng = val.Get ();
gjc@4249:         NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2:0.15");
gjc@4249:       }
mathieu@4243:     }
mathieu@4243: 
mathieu@4243:     return result;
mathieu@4243:   }
mathieu@4243: };
mathieu@4243: 
mathieu@4243: 
mathieu@4243: static RandomVariableTest g_random_variable_tests;
mathieu@4243: 
mathieu@4243: }//namespace ns3
mathieu@4243: 
mathieu@4243: #endif /* RUN_SELF_TESTS */