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