1 /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */

     2 //

     3 // Copyright (c) 2006 Georgia Tech Research Corporation

     4 //

     5 // This program is free software; you can redistribute it and/or modify

     6 // it under the terms of the GNU General Public License version 2 as

     7 // published by the Free Software Foundation;

     8 //

     9 // This program is distributed in the hope that it will be useful,

    10 // but WITHOUT ANY WARRANTY; without even the implied warranty of

    11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

    12 // GNU General Public License for more details.

    13 //

    14 // You should have received a copy of the GNU General Public License

    15 // along with this program; if not, write to the Free Software

    16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    17 //

    18 // Author: Rajib Bhattacharjea<raj.b@gatech.edu>

    19 //

    20 

    21 #include <iostream>

    22 

    23 #include <math.h>

    24 #include <stdlib.h>

    25 #include <sys/time.h>			// for gettimeofday

    26 #include <unistd.h>

    27 #include <iostream>

    28 #include <sys/types.h>

    29 #include <sys/stat.h>

    30 #include <fcntl.h>       

    31 

    32 

    33 #include "random-variable.h"

    34 #include "rng-stream.h"

    35 #include "fatal-error.h"

    36 

    37 using namespace std;

    38 

    39 namespace ns3{

    40 // Seed methods

    41 

    42 Seed::~Seed()

    43 {

    44 }

    45 

    46 RandomSeed::RandomSeed()

    47 {

    48 }

    49 

    50 RandomSeed::~RandomSeed()

    51 {

    52 }

    53 

    54 bool RandomSeed::IsRandom() const 

    55 {

    56   return true;

    57 }

    58 

    59 ConstantSeed::~ConstantSeed()

    60 {

    61 }

    62 

    63 bool ConstantSeed::IsRandom() const 

    64 {

    65   return false;

    66 }

    67 

    68 ConstantSeed::ConstantSeed(uint32_t s)

    69 {

    70   seeds[0] = s;

    71   seeds[1] = s;

    72   seeds[2] = s;

    73   seeds[3] = s;

    74   seeds[4] = s;

    75   seeds[5] = s;

    76 }

    77 

    78 ConstantSeed::ConstantSeed(uint32_t s0, uint32_t s1, uint32_t s2,

    79                            uint32_t s3, uint32_t s4, uint32_t s5)

    80 {

    81   seeds[0] = s0;

    82   seeds[1] = s1;

    83   seeds[2] = s2;

    84   seeds[3] = s3;

    85   seeds[4] = s4;

    86   seeds[5] = s5;

    87 }

    88 //-----------------------------------------------------------------------------

    89 //-----------------------------------------------------------------------------

    90 // RandomVariable methods

    91 

    92 bool          RandomVariable::initialized = false;   // True if RngStream seed set 

    93 bool          RandomVariable::useDevRandom = false;  // True if use /dev/random desired

    94 bool          RandomVariable::globalSeedSet = false; // True if GlobalSeed called

    95 int           RandomVariable::devRandom = -1;

    96 uint32_t        RandomVariable::globalSeed[6];

    97 unsigned long RandomVariable::heuristic_sequence;

    98 

    99 RandomVariable::RandomVariable() 

   100 {

   101   m_generator = new RngStream();

   102   RandomVariable::Initialize(); // sets the seed for the static object

   103   m_generator->InitializeStream();

   104 }

   105 

   106 RandomVariable::~RandomVariable()

   107 {

   108   delete m_generator;

   109 }

   110 

   111 uint32_t RandomVariable::GetIntValue() 

   112 {

   113   return (uint32_t)GetValue();

   114 }

   115 

   116 void RandomVariable::UseDevRandom(bool udr) 

   117 {

   118   RandomVariable::useDevRandom = udr;

   119 }

   120 

   121 bool RandomVariable::SetSeed(const Seed& s)

   122 {

   123   // Seed this stream with the specified seed

   124   if (s.IsRandom())

   125     {

   126       uint32_t seeds[6];

   127       while(true)

   128         { // Insure seeds are valid

   129           GetRandomSeeds(seeds);

   130           if (RngStream::CheckSeed(seeds)) break;

   131         }

   132       m_generator->SetSeeds(seeds);

   133       return true;

   134     }

   135   // Not random seed, use specified

   136   const ConstantSeed& cs = (ConstantSeed&)s;

   137   if (!RngStream::CheckSeed(cs.seeds))

   138     {

   139       cout << "Constant seed failed valid check" << endl;

   140       return false; // Seed is not valid

   141     }

   142   m_generator->SetSeeds(cs.seeds);

   143   return true;

   144 }

   145 

   146 //-----------------------------------------------------------------------------

   147 //-----------------------------------------------------------------------------

   148 // RandomVariable static methods

   149 void RandomVariable::UseGlobalSeed(const Seed& s)

   150 {

   151   if (RandomVariable::globalSeedSet)

   152     {

   153       cout << "Random number generator already initialized!" << endl;

   154       cout << "Call to RandomVariable::UseGlobalSeed() ignored" << endl;

   155       return;

   156     }

   157   if (s.IsRandom()) return; // Random seed is the default

   158   const ConstantSeed& cs = (ConstantSeed&)s;

   159   RandomVariable::globalSeed[0] = cs.seeds[0];

   160   RandomVariable::globalSeed[1] = cs.seeds[1];

   161   RandomVariable::globalSeed[2] = cs.seeds[2];

   162   RandomVariable::globalSeed[3] = cs.seeds[3];

   163   RandomVariable::globalSeed[4] = cs.seeds[4];

   164   RandomVariable::globalSeed[5] = cs.seeds[5];

   165   if (!RngStream::CheckSeed(RandomVariable::globalSeed))

   166   	NS_FATAL_ERROR("Invalid seed");

   167   

   168   RandomVariable::globalSeedSet = true;

   169 }

   170 

   171 void RandomVariable::Initialize()

   172 { 

   173   if (RandomVariable::initialized) return; // Already initialized and seeded

   174   RandomVariable::initialized = true;

   175   if (!RandomVariable::globalSeedSet)

   176     { // No global seed, try a random one

   177       GetRandomSeeds(globalSeed);

   178     }

   179   // Seed the RngStream package

   180   RngStream::SetPackageSeed(globalSeed);

   181 }

   182 

   183 void RandomVariable::GetRandomSeeds(uint32_t seeds[6])

   184 {

   185   // Check if /dev/random exists

   186   if (RandomVariable::useDevRandom && RandomVariable::devRandom < 0)

   187     {

   188       RandomVariable::devRandom = open("/dev/random", O_RDONLY);

   189     }

   190   if (RandomVariable::devRandom > 0)

   191     { // Use /dev/random

   192       while(true)

   193         {

   194           for (int i = 0; i < 6; ++i)

   195             {

   196               read(RandomVariable::devRandom, &seeds[i], sizeof(seeds[i]));

   197             }

   198           if (RngStream::CheckSeed(seeds)) break; // Got a valid one

   199         }

   200     }

   201   else

   202     { // Seed from time of day (code borrowed from ns2 random seeding)

   203       // Thanks to John Heidemann for this technique

   204       while(true)

   205         {

   206           timeval tv;

   207           gettimeofday(&tv, 0);

   208           seeds[0] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   209               & 0x7fffffff;

   210           gettimeofday(&tv, 0);

   211           seeds[1] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   212               & 0x7fffffff;

   213           gettimeofday(&tv, 0);

   214           seeds[2] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   215               & 0x7fffffff;

   216           gettimeofday(&tv, 0);

   217           seeds[3] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   218               & 0x7fffffff;

   219           gettimeofday(&tv, 0);

   220           seeds[4] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   221               & 0x7fffffff;

   222           gettimeofday(&tv, 0);

   223           seeds[5] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))

   224               & 0x7fffffff;

   225           if (RngStream::CheckSeed(seeds)) break; // Got a valid one

   226         }

   227     }

   228 }

   229 

   230 //-----------------------------------------------------------------------------

   231 //-----------------------------------------------------------------------------

   232 // UniformVariable methods

   233 UniformVariable::UniformVariable() 

   234   : m_min(0), m_max(1.0) { }

   235   

   236 UniformVariable::UniformVariable(double s, double l) 

   237   : m_min(s), m_max(l) { }

   238 

   239 UniformVariable::UniformVariable(const UniformVariable& c) 

   240   : m_min(c.m_min), m_max(c.m_max) { }

   241 

   242 double UniformVariable::GetValue()

   243 {

   244   return m_min + m_generator->RandU01() * (m_max - m_min);

   245 }

   246 

   247 RandomVariable* UniformVariable::Copy() const

   248 {

   249   return new UniformVariable(*this);

   250 }

   251 //-----------------------------------------------------------------------------

   252 //-----------------------------------------------------------------------------

   253 // ConstantVariable methods

   254 ConstantVariable::ConstantVariable() 

   255   : m_const(0) { }

   256 

   257 ConstantVariable::ConstantVariable(double c) 

   258   : m_const(c) { };

   259   

   260 ConstantVariable::ConstantVariable(const ConstantVariable& c) 

   261   : m_const(c.m_const) { }

   262 

   263 void ConstantVariable::NewConstant(double c) 

   264   { m_const = c;}

   265   

   266 double ConstantVariable::GetValue()

   267 {

   268   return m_const;

   269 }

   270 

   271 uint32_t ConstantVariable::GetIntValue()

   272 {

   273   return (uint32_t)m_const;

   274 }

   275 

   276 RandomVariable* ConstantVariable::Copy() const

   277 {

   278   return new ConstantVariable(*this);

   279 }

   280 //-----------------------------------------------------------------------------

   281 //-----------------------------------------------------------------------------

   282 // SequentialVariable methods

   283 SequentialVariable::SequentialVariable(double f, double l, double i, uint32_t c)

   284   : m_min(f), m_max(l), m_increment(ConstantVariable(i).Copy()), m_consecutive(c),

   285     m_current(f), m_currentConsecutive(0)

   286 {

   287 }

   288 

   289 SequentialVariable::SequentialVariable(double f, double l, const RandomVariable& i, uint32_t c)

   290   : m_min(f), m_max(l), m_increment(i.Copy()), m_consecutive(c),

   291     m_current(f), m_currentConsecutive(0)

   292 {

   293 }

   294 

   295 SequentialVariable::SequentialVariable(const SequentialVariable& c)

   296   : m_min(c.m_min), m_max(c.m_max),

   297     m_increment(c.m_increment->Copy()), m_consecutive(c.m_consecutive),

   298     m_current(c.m_current), m_currentConsecutive(c.m_currentConsecutive)

   299 {

   300 }

   301 

   302 double SequentialVariable::GetValue()

   303 { // Return a sequential series of values

   304   double r = m_current;

   305   if (++m_currentConsecutive == m_consecutive)

   306     { // Time to advance to next

   307       m_currentConsecutive = 0;

   308       m_current += m_increment->GetValue();

   309       if (m_current >= m_max)

   310         m_current = m_min + (m_current - m_max);

   311     }

   312   return r;

   313 }

   314 

   315 RandomVariable* SequentialVariable::Copy() const

   316 {

   317   return new SequentialVariable(*this);

   318 }

   319 //-----------------------------------------------------------------------------

   320 //-----------------------------------------------------------------------------

   321 // ExponentialVariable methods

   322 ExponentialVariable::ExponentialVariable() 

   323   : m_mean(1.0), m_bound(0) { }

   324   

   325 ExponentialVariable::ExponentialVariable(double m) 

   326   : m_mean(m), m_bound(0) { }

   327   

   328 ExponentialVariable::ExponentialVariable(double m, double b) 

   329   : m_mean(m), m_bound(b) { }

   330   

   331 ExponentialVariable::ExponentialVariable(const ExponentialVariable& c) 

   332   : m_mean(c.m_mean), m_bound(c.m_bound) { }

   333 

   334 double ExponentialVariable::GetValue()

   335 {

   336   double r = -m_mean*log(m_generator->RandU01());

   337   if (m_bound != 0 && r > m_bound) return m_bound;

   338   return r;

   339 }

   340 

   341 RandomVariable* ExponentialVariable::Copy() const

   342 {

   343   return new ExponentialVariable(*this);

   344 }

   345 //-----------------------------------------------------------------------------

   346 //-----------------------------------------------------------------------------

   347 // ParetoVariable methods

   348 ParetoVariable::ParetoVariable() 

   349   : m_mean(1.0), m_shape(1.5), m_bound(0) { }

   350 

   351 ParetoVariable::ParetoVariable(double m) 

   352   : m_mean(m), m_shape(1.5), m_bound(0) { }

   353 

   354 ParetoVariable::ParetoVariable(double m, double s) 

   355     : m_mean(m), m_shape(s), m_bound(0) { }

   356 

   357 ParetoVariable::ParetoVariable(double m, double s, double b) 

   358   : m_mean(m), m_shape(s), m_bound(b) { }

   359 

   360 ParetoVariable::ParetoVariable(const ParetoVariable& c) 

   361   : m_mean(c.m_mean), m_shape(c.m_shape), m_bound(c.m_bound) { }

   362 

   363 double ParetoVariable::GetValue()

   364 {

   365   double scale = m_mean * ( m_shape - 1.0) / m_shape;

   366   double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));

   367   if (m_bound != 0 && r > m_bound) return m_bound;

   368   return r;

   369 }

   370 

   371 RandomVariable* ParetoVariable::Copy() const

   372 {

   373   return new ParetoVariable(*this);

   374 }

   375 //-----------------------------------------------------------------------------

   376 //-----------------------------------------------------------------------------

   377 // WeibullVariable methods

   378 WeibullVariable::WeibullVariable() : m_mean(1.0), m_alpha(1), m_bound(0) { }

   379 WeibullVariable::WeibullVariable(double m) 

   380   : m_mean(m), m_alpha(1), m_bound(0) { }

   381 WeibullVariable::WeibullVariable(double m, double s) 

   382   : m_mean(m), m_alpha(s), m_bound(0) { }

   383 WeibullVariable::WeibullVariable(double m, double s, double b) 

   384   : m_mean(m), m_alpha(s), m_bound(b) { };

   385 WeibullVariable::WeibullVariable(const WeibullVariable& c) 

   386   : m_mean(c.m_mean), m_alpha(c.m_alpha), m_bound(c.m_bound) { }

   387 

   388 double WeibullVariable::GetValue()

   389 {

   390   double exponent = 1.0 / m_alpha;

   391   double r = m_mean * pow( -log(m_generator->RandU01()), exponent);

   392   if (m_bound != 0 && r > m_bound) return m_bound;

   393   return r;

   394 }

   395 

   396 RandomVariable* WeibullVariable::Copy() const

   397 {

   398   return new WeibullVariable(*this);

   399 }

   400 //-----------------------------------------------------------------------------

   401 //-----------------------------------------------------------------------------

   402 // NormalVariable methods

   403 NormalVariable::NormalVariable() 

   404   : m_mean(0.0), m_variance(1.0), m_bound(INFINITE_VALUE), m_nextValid(false){}

   405 

   406 NormalVariable::NormalVariable(double m, double v, double b)

   407   : m_mean(m), m_variance(v), m_bound(b), m_nextValid(false) { }

   408 

   409 NormalVariable::NormalVariable(const NormalVariable& c)

   410   : m_mean(c.m_mean), m_variance(c.m_variance), m_bound(c.m_bound) { }

   411 

   412 double NormalVariable::GetValue()

   413 {

   414   if (m_nextValid)

   415     { // use previously generated

   416       m_nextValid = false;

   417       return m_next;

   418     }

   419   while(1)

   420     { // See Simulation Modeling and Analysis p. 466 (Averill Law)

   421       // for algorithm

   422       double u1 = m_generator->RandU01();

   423       double u2 = m_generator->RandU01();;

   424       double v1 = 2 * u1 - 1;

   425       double v2 = 2 * u2 - 1;

   426       double w = v1 * v1 + v2 * v2;

   427       if (w <= 1.0)

   428         { // Got good pair

   429           double y = sqrt((-2 * log(w))/w);

   430           m_next = m_mean + v2 * y * sqrt(m_variance);

   431           if (fabs(m_next) > m_bound) m_next = m_bound * (m_next)/fabs(m_next);

   432           m_nextValid = true;

   433           double x1 = m_mean + v1 * y * sqrt(m_variance);

   434           if (fabs(x1) > m_bound) x1 = m_bound * (x1)/fabs(x1);

   435           return x1;

   436         }

   437     }

   438 }

   439 

   440 RandomVariable* NormalVariable::Copy() const

   441 {

   442   return new NormalVariable(*this);

   443 }

   444 

   445 //-----------------------------------------------------------------------------

   446 //-----------------------------------------------------------------------------

   447 // ValueCDF methods

   448 ValueCDF::ValueCDF() 

   449   : value(0.0), cdf(0.0){ }

   450 ValueCDF::ValueCDF(double v, double c) 

   451   : value(v), cdf(c) { }

   452 ValueCDF::ValueCDF(const ValueCDF& c) 

   453   : value(c.value), cdf(c.cdf) { }

   454 

   455 //-----------------------------------------------------------------------------

   456 //-----------------------------------------------------------------------------

   457 // EmpiricalVariable methods

   458 EmpiricalVariable::EmpiricalVariable() 

   459   : validated(false) { }

   460 

   461 EmpiricalVariable::EmpiricalVariable(const EmpiricalVariable& c) 

   462   : validated(c.validated), emp(c.emp) { }

   463 

   464 EmpiricalVariable::~EmpiricalVariable() { }

   465 

   466 double EmpiricalVariable::GetValue()

   467 { // Return a value from the empirical distribution

   468   // This code based (loosely) on code by Bruce Mah (Thanks Bruce!)

   469   if (emp.size() == 0) return 0.0; // HuH? No empirical data

   470   if (!validated) Validate();      // Insure in non-decreasing

   471   double r = m_generator->RandU01();

   472   if (r <= emp.front().cdf)return emp.front().value; // Less than first

   473   if (r >= emp.back().cdf) return emp.back().value;  // Greater than last

   474   // Binary search

   475   std::vector<ValueCDF>::size_type bottom = 0;

   476   std::vector<ValueCDF>::size_type top = emp.size() - 1;

   477   while(1)

   478     {

   479       std::vector<ValueCDF>::size_type c = (top + bottom) / 2;

   480       if (r >= emp[c].cdf && r < emp[c+1].cdf)

   481         { // Found it

   482           return Interpolate(emp[c].cdf, emp[c+1].cdf,

   483                              emp[c].value, emp[c+1].value,

   484                              r);

   485         }

   486       // Not here, adjust bounds

   487       if (r < emp[c].cdf) top    = c - 1;

   488       else                bottom = c + 1;

   489     }

   490 }

   491 

   492 RandomVariable* EmpiricalVariable::Copy() const

   493 {

   494   return new EmpiricalVariable(*this);

   495 }

   496 

   497 void EmpiricalVariable::CDF(double v, double c)

   498 { // Add a new empirical datapoint to the empirical cdf

   499   // NOTE.   These MUST be inserted in non-decreasing order

   500   emp.push_back(ValueCDF(v, c));

   501 }

   502 

   503 void EmpiricalVariable::Validate()

   504 {

   505   ValueCDF prior;

   506   for (std::vector<ValueCDF>::size_type i = 0; i < emp.size(); ++i)

   507     {

   508       ValueCDF& current = emp[i];

   509       if (current.value < prior.value || current.cdf < prior.cdf)

   510         { // Error

   511           cout << "Empirical Dist error,"

   512                << " current value " << current.value

   513                << " prior value "   << prior.value

   514                << " current cdf "   << current.cdf

   515                << " prior cdf "     << prior.cdf << endl;

   516           NS_FATAL_ERROR("Empirical Dist error");

   517         }

   518       prior = current;

   519     }

   520   validated = true;

   521 }

   522 

   523 double EmpiricalVariable::Interpolate(double c1, double c2,

   524                                 double v1, double v2, double r)

   525 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)

   526   return (v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));

   527 }

   528 

   529 //-----------------------------------------------------------------------------

   530 //-----------------------------------------------------------------------------

   531 // Integer EmpiricalVariable methods

   532 IntEmpiricalVariable::IntEmpiricalVariable() { }

   533 

   534 uint32_t IntEmpiricalVariable::GetIntValue()

   535 {

   536   return (uint32_t)GetValue();

   537 }

   538 

   539 RandomVariable* IntEmpiricalVariable::Copy() const

   540 {

   541   return new IntEmpiricalVariable(*this);

   542 }

   543 

   544 

   545 double IntEmpiricalVariable::Interpolate(double c1, double c2,

   546                                    double v1, double v2, double r)

   547 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)

   548   return ceil(v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));

   549 }

   550 

   551 

   552 //-----------------------------------------------------------------------------

   553 //-----------------------------------------------------------------------------

   554 // DeterministicVariable

   555 DeterministicVariable::DeterministicVariable(double* d, uint32_t c)

   556     : count(c), next(c), data(d)

   557 { // Nothing else needed

   558 }

   559 

   560 DeterministicVariable::~DeterministicVariable() { }

   561   

   562 double DeterministicVariable::GetValue()

   563 {

   564   if (next == count) next = 0;

   565   return data[next++];

   566 }

   567 

   568 RandomVariable* DeterministicVariable::Copy() const

   569 {

   570   return new DeterministicVariable(*this);

   571 }

   572 

   573 }//namespace ns3

   574