View | Details | Raw Unified | Return to bug 578
Collapse All | Expand All

(-)a/src/core/random-variable.cc (-552 / +333 lines)
 Lines 16-21    Link Here 
16
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
16
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17
//
17
//
18
// Author: Rajib Bhattacharjea<raj.b@gatech.edu>
18
// Author: Rajib Bhattacharjea<raj.b@gatech.edu>
19
// Author: Hadi Arbabi<marbabi@cs.odu.edu>
19
//
20
//
20
21
21
#include <iostream>
22
#include <iostream>
 Lines 32-37    Link Here 
32
33
33
34
34
#include "assert.h"
35
#include "assert.h"
36
#include "config.h"
37
#include "integer.h"
35
#include "random-variable.h"
38
#include "random-variable.h"
36
#include "rng-stream.h"
39
#include "rng-stream.h"
37
#include "fatal-error.h"
40
#include "fatal-error.h"
 Lines 41-46    Link Here 
41
namespace ns3{
44
namespace ns3{
42
45
43
//-----------------------------------------------------------------------------
46
//-----------------------------------------------------------------------------
47
// Seed Manager
48
//-----------------------------------------------------------------------------
49
50
uint32_t SeedManager::GetSeed()
51
{
52
  uint32_t s[6];
53
  RngStream::GetPackageSeed (s);
54
  NS_ASSERT(
55
              s[0] == s[1] &&
56
              s[0] == s[2] &&
57
              s[0] == s[3] &&
58
              s[0] == s[4] &&
59
              s[0] == s[5]    
60
            );
61
  return s[0];
62
}
63
64
void SeedManager::SetSeed(uint32_t seed)
65
{
66
  Config::SetGlobal("RngSeed", IntegerValue(seed));
67
}
68
69
void SeedManager::SetRun(uint32_t run)
70
{
71
  Config::SetGlobal("RngRun", IntegerValue(run));
72
}
73
74
uint32_t SeedManager::GetRun()
75
{
76
  return RngStream::GetPackageRun ();
77
}
78
79
bool SeedManager::CheckSeed (uint32_t seed)
80
{
81
  return RngStream::CheckSeed(seed);
82
}
83
84
//-----------------------------------------------------------------------------
44
//-----------------------------------------------------------------------------
85
//-----------------------------------------------------------------------------
45
// RandomVariableBase methods
86
// RandomVariableBase methods
46
87
 Lines 54-125    Link Here 
54
  virtual double  GetValue() = 0;
95
  virtual double  GetValue() = 0;
55
  virtual uint32_t GetInteger();
96
  virtual uint32_t GetInteger();
56
  virtual RandomVariableBase*   Copy(void) const = 0;
97
  virtual RandomVariableBase*   Copy(void) const = 0;
57
  virtual void GetSeed(uint32_t seed[6]);
58
98
59
  static  void UseDevRandom(bool udr = true);
60
  static void UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, 
61
                            uint32_t s3, uint32_t s4, uint32_t s5);
62
  static void SetRunNumber(uint32_t n);
63
private:
64
  static void GetRandomSeeds(uint32_t seeds[6]);
65
private:
66
  static bool useDevRandom;    // True if using /dev/random desired
67
  static bool globalSeedSet;   // True if global seed has been specified
68
  static int  devRandom;       // File handle for /dev/random
69
  static uint32_t globalSeed[6]; // The global seed to use
70
  friend class RandomVariableInitializer;
71
protected:
99
protected:
72
  static unsigned long heuristic_sequence;
73
  static RngStream* m_static_generator;
74
  static uint32_t runNumber;
75
  static void Initialize();    // Initialize  the RNG system
76
  static bool initialized;     // True if package seed is set 
77
  RngStream* m_generator;  //underlying generator being wrapped
100
  RngStream* m_generator;  //underlying generator being wrapped
78
};
101
};
79
102
80
81
82
bool          RandomVariableBase::initialized = false;   // True if RngStream seed set 
83
bool          RandomVariableBase::useDevRandom = false;  // True if use /dev/random
84
bool          RandomVariableBase::globalSeedSet = false; // True if GlobalSeed called
85
int           RandomVariableBase::devRandom = -1;
86
uint32_t      RandomVariableBase::globalSeed[6];
87
unsigned long RandomVariableBase::heuristic_sequence;
88
RngStream*    RandomVariableBase::m_static_generator = 0;
89
uint32_t      RandomVariableBase::runNumber = 0;
90
91
//the static object random_variable_initializer initializes the static members
92
//of RandomVariable
93
static class RandomVariableInitializer
94
{
95
  public:
96
  RandomVariableInitializer()
97
  {
98
//     RandomVariableBase::Initialize(); // sets the static package seed
99
//     RandomVariableBase::m_static_generator = new RngStream();
100
//     RandomVariableBase::m_static_generator->InitializeStream();
101
  }
102
  ~RandomVariableInitializer()
103
  {
104
    delete RandomVariableBase::m_static_generator;
105
  }
106
} random_variable_initializer;
107
108
RandomVariableBase::RandomVariableBase() 
103
RandomVariableBase::RandomVariableBase() 
109
  : m_generator(NULL)
104
  : m_generator(NULL)
110
{
105
{
111
//   m_generator = new RngStream();
112
//   m_generator->InitializeStream();
113
//   m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
114
}
106
}
115
107
116
RandomVariableBase::RandomVariableBase(const RandomVariableBase& r)
108
RandomVariableBase::RandomVariableBase(const RandomVariableBase& r)
117
  :m_generator(0)
109
  :m_generator(0)
118
{
110
{
119
  if(r.m_generator)
111
  if (r.m_generator)
120
  {
112
    {
121
    m_generator = new RngStream(*r.m_generator);
113
      m_generator = new RngStream(*r.m_generator);
122
  }
114
    }
123
}
115
}
124
116
125
RandomVariableBase::~RandomVariableBase()
117
RandomVariableBase::~RandomVariableBase()
 Lines 132-246    Link Here 
132
  return (uint32_t)GetValue();
124
  return (uint32_t)GetValue();
133
}
125
}
134
126
135
void RandomVariableBase::UseDevRandom(bool udr) 
127
//-------------------------------------------------------
136
{
137
  RandomVariableBase::useDevRandom = udr;
138
}
139
140
void RandomVariableBase::GetSeed(uint32_t seed[6])
141
{
142
  if(!m_generator)
143
  {
144
    m_generator = new RngStream();
145
    m_generator->InitializeStream();
146
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
147
  }
148
  m_generator->GetState(seed);
149
}
150
151
//-----------------------------------------------------------------------------
152
//-----------------------------------------------------------------------------
153
// RandomVariableBase static methods
154
void RandomVariableBase::UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, 
155
                                   uint32_t s3, uint32_t s4, uint32_t s5)
156
{
157
  if (RandomVariableBase::globalSeedSet)
158
    {
159
      NS_FATAL_ERROR ("Random number generator already initialized! "
160
                      "Call to RandomVariableBase::UseGlobalSeed() ignored");
161
    }
162
  RandomVariableBase::globalSeed[0] = s0;
163
  RandomVariableBase::globalSeed[1] = s1;
164
  RandomVariableBase::globalSeed[2] = s2;
165
  RandomVariableBase::globalSeed[3] = s3;
166
  RandomVariableBase::globalSeed[4] = s4;
167
  RandomVariableBase::globalSeed[5] = s5;
168
  if (!RngStream::CheckSeed(RandomVariableBase::globalSeed))
169
    NS_FATAL_ERROR("Invalid seed");
170
  
171
  RandomVariableBase::globalSeedSet = true;
172
}
173
174
void RandomVariableBase::Initialize()
175
{ 
176
  if (RandomVariableBase::initialized) return; // Already initialized and seeded
177
  RandomVariableBase::initialized = true;
178
  if (!RandomVariableBase::globalSeedSet)
179
    { // No global seed, try a random one
180
      GetRandomSeeds(globalSeed);
181
    }
182
  // Seed the RngStream package
183
  RngStream::SetPackageSeed(globalSeed);
184
}
185
186
void RandomVariableBase::GetRandomSeeds(uint32_t seeds[6])
187
{
188
  // Check if /dev/random exists
189
  if (RandomVariableBase::useDevRandom && RandomVariableBase::devRandom < 0)
190
    {
191
      RandomVariableBase::devRandom = open("/dev/random", O_RDONLY);
192
    }
193
  if (RandomVariableBase::devRandom > 0)
194
    { // Use /dev/random
195
      while(true)
196
        {
197
          for (int i = 0; i < 6; ++i)
198
            {
199
              ssize_t bytes_read = read (RandomVariableBase::devRandom,
200
                                         &seeds[i], sizeof (seeds[i]));
201
              if (bytes_read != sizeof (seeds[i]))
202
                {
203
                  NS_FATAL_ERROR ("Read from /dev/random failed");
204
                }
205
            }
206
          if (RngStream::CheckSeed(seeds)) break; // Got a valid one
207
        }
208
    }
209
  else
210
    { // Seed from time of day (code borrowed from ns2 random seeding)
211
      // Thanks to John Heidemann for this technique
212
      while(true)
213
        {
214
          timeval tv;
215
          gettimeofday(&tv, 0);
216
          seeds[0] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
217
              & 0x7fffffff;
218
          gettimeofday(&tv, 0);
219
          seeds[1] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
220
              & 0x7fffffff;
221
          gettimeofday(&tv, 0);
222
          seeds[2] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
223
              & 0x7fffffff;
224
          gettimeofday(&tv, 0);
225
          seeds[3] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
226
              & 0x7fffffff;
227
          gettimeofday(&tv, 0);
228
          seeds[4] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
229
              & 0x7fffffff;
230
          gettimeofday(&tv, 0);
231
          seeds[5] = (tv.tv_sec^tv.tv_usec^(++heuristic_sequence <<8))
232
              & 0x7fffffff;
233
          if (RngStream::CheckSeed(seeds)) break; // Got a valid one
234
        }
235
    }
236
}
237
238
void RandomVariableBase::SetRunNumber(uint32_t n)
239
{
240
  runNumber = n;
241
}
242
243
244
128
245
RandomVariable::RandomVariable()
129
RandomVariable::RandomVariable()
246
  : m_variable (0)
130
  : m_variable (0)
 Lines 277-309    Link Here 
277
{
161
{
278
  return m_variable->GetInteger ();
162
  return m_variable->GetInteger ();
279
}
163
}
280
void 
164
281
RandomVariable::GetSeed(uint32_t seed[6]) const
282
{
283
  return m_variable->GetSeed (seed);
284
}
285
void 
286
RandomVariable::UseDevRandom(bool udr)
287
{
288
  RandomVariableBase::UseDevRandom (udr);
289
}
290
void 
291
RandomVariable::UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, 
292
                              uint32_t s3, uint32_t s4, uint32_t s5)
293
{
294
  RandomVariableBase::UseGlobalSeed (s0, s1, s2, s3, s4, s5);
295
}
296
void 
297
RandomVariable::SetRunNumber(uint32_t n)
298
{
299
  RandomVariableBase::SetRunNumber (n);
300
}
301
RandomVariableBase *
165
RandomVariableBase *
302
RandomVariable::Peek (void) const
166
RandomVariable::Peek (void) const
303
{
167
{
304
  return m_variable;
168
  return m_variable;
305
}
169
}
306
170
171
307
ATTRIBUTE_VALUE_IMPLEMENT (RandomVariable);
172
ATTRIBUTE_VALUE_IMPLEMENT (RandomVariable);
308
ATTRIBUTE_CHECKER_IMPLEMENT (RandomVariable);
173
ATTRIBUTE_CHECKER_IMPLEMENT (RandomVariable);
309
174
 Lines 335-349    Link Here 
335
   * \return A value between low and high values specified by the constructor
200
   * \return A value between low and high values specified by the constructor
336
   */
201
   */
337
  virtual double GetValue();
202
  virtual double GetValue();
203
204
  /**
205
   * \return A value between low and high values specified by parameters
206
   */
207
  virtual double GetValue(double s, double l);
208
  
338
  virtual RandomVariableBase*  Copy(void) const;
209
  virtual RandomVariableBase*  Copy(void) const;
339
210
340
public:
341
  /**
342
   * \param s Low end of the range
343
   * \param l High end of the range
344
   * \return A uniformly distributed random number between s and l
345
   */
346
  static double GetSingleValue(double s, double l);
347
private:
211
private:
348
  double m_min;
212
  double m_min;
349
  double m_max;
213
  double m_max;
 Lines 372-418    Link Here 
372
236
373
double UniformVariableImpl::GetValue()
237
double UniformVariableImpl::GetValue()
374
{
238
{
375
  if(!RandomVariableBase::initialized)
376
  {
377
    RandomVariableBase::Initialize();
378
  }
379
  if(!m_generator)
239
  if(!m_generator)
380
  {
240
    {
381
    m_generator = new RngStream();
241
      m_generator = new RngStream();
382
    m_generator->InitializeStream();
242
    }
383
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
384
  }
385
  return m_min + m_generator->RandU01() * (m_max - m_min);
243
  return m_min + m_generator->RandU01() * (m_max - m_min);
386
}
244
}
387
245
246
double UniformVariableImpl::GetValue(double s, double l) 
247
{
248
  if(!m_generator)
249
    {
250
      m_generator = new RngStream();
251
    }
252
    return s + m_generator->RandU01() * (l-s); 
253
}
254
388
RandomVariableBase* UniformVariableImpl::Copy() const
255
RandomVariableBase* UniformVariableImpl::Copy() const
389
{
256
{
390
  return new UniformVariableImpl(*this);
257
  return new UniformVariableImpl(*this);
391
}
258
}
392
259
393
double UniformVariableImpl::GetSingleValue(double s, double l)
394
{
395
  if(!RandomVariableBase::m_static_generator)
396
  {
397
    RandomVariableBase::Initialize(); // sets the static package seed
398
    RandomVariableBase::m_static_generator = new RngStream();
399
    RandomVariableBase::m_static_generator->InitializeStream();
400
  }
401
  return s + m_static_generator->RandU01() * (l - s);;
402
}
403
404
UniformVariable::UniformVariable()
260
UniformVariable::UniformVariable()
405
  : RandomVariable (UniformVariableImpl ())
261
  : RandomVariable (UniformVariableImpl ())
406
{}
262
{}
407
UniformVariable::UniformVariable(double s, double l)
263
UniformVariable::UniformVariable(double s, double l)
408
  : RandomVariable (UniformVariableImpl (s, l))
264
  : RandomVariable (UniformVariableImpl (s, l))
409
{}
265
{}
410
double 
266
411
UniformVariable::GetSingleValue(double s, double l)
267
double UniformVariable::GetValue(void) const
412
{
268
{
413
  return UniformVariableImpl::GetSingleValue (s, l);
269
  return this->RandomVariable::GetValue ();
414
}
270
}
415
271
272
double UniformVariable::GetValue(double s, double l)
273
{
274
  return ((UniformVariableImpl*)Peek())->GetValue(s,l);
275
}
276
277
uint32_t UniformVariable::GetInteger (uint32_t s, uint32_t l)
278
{
279
  NS_ASSERT(s <= l);
280
  return static_cast<uint32_t>( GetValue(s, l+1) );
281
}
416
282
417
//-----------------------------------------------------------------------------
283
//-----------------------------------------------------------------------------
418
//-----------------------------------------------------------------------------
284
//-----------------------------------------------------------------------------
 Lines 623-635    Link Here 
623
   */
489
   */
624
  virtual double GetValue();
490
  virtual double GetValue();
625
  virtual RandomVariableBase* Copy(void) const;
491
  virtual RandomVariableBase* Copy(void) const;
626
public:
492
627
  /**
628
   * \param m The mean of the distribution from which the return value is drawn
629
   * \param b The upper bound value desired, beyond which values get clipped
630
   * \return A random number from an exponential distribution with mean m
631
   */
632
  static double GetSingleValue(double m, double b=0);
633
private:
493
private:
634
  double m_mean;  // Mean value of RV
494
  double m_mean;  // Mean value of RV
635
  double m_bound; // Upper bound on value (if non-zero)
495
  double m_bound; // Upper bound on value (if non-zero)
 Lines 649-685    Link Here 
649
509
650
double ExponentialVariableImpl::GetValue()
510
double ExponentialVariableImpl::GetValue()
651
{
511
{
652
  if(!RandomVariableBase::initialized)
653
  {
654
    RandomVariableBase::Initialize();
655
  }
656
  if(!m_generator)
512
  if(!m_generator)
657
  {
513
    {
658
    m_generator = new RngStream();
514
      m_generator = new RngStream();
659
    m_generator->InitializeStream();
515
    }
660
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
516
  while(1)
661
  }
517
    {
662
  double r = -m_mean*log(m_generator->RandU01());
518
      double r = -m_mean*log(m_generator->RandU01());
663
  if (m_bound != 0 && r > m_bound) return m_bound;
519
      if (m_bound == 0 || r <= m_bound) return r;
664
  return r;
520
      //otherwise, try again
521
    }
665
}
522
}
666
523
667
RandomVariableBase* ExponentialVariableImpl::Copy() const
524
RandomVariableBase* ExponentialVariableImpl::Copy() const
668
{
525
{
669
  return new ExponentialVariableImpl(*this);
526
  return new ExponentialVariableImpl(*this);
670
}
527
}
671
double ExponentialVariableImpl::GetSingleValue(double m, double b/*=0*/)
672
{
673
  if(!RandomVariableBase::m_static_generator)
674
  {
675
    RandomVariableBase::Initialize(); // sets the static package seed
676
    RandomVariableBase::m_static_generator = new RngStream();
677
    RandomVariableBase::m_static_generator->InitializeStream();
678
  }
679
  double r = -m*log(m_static_generator->RandU01());
680
  if (b != 0 && r > b) return b;
681
  return r;
682
}
683
528
684
ExponentialVariable::ExponentialVariable()
529
ExponentialVariable::ExponentialVariable()
685
  : RandomVariable (ExponentialVariableImpl ())
530
  : RandomVariable (ExponentialVariableImpl ())
 Lines 690-700    Link Here 
690
ExponentialVariable::ExponentialVariable(double m, double b)
535
ExponentialVariable::ExponentialVariable(double m, double b)
691
  : RandomVariable (ExponentialVariableImpl (m, b))
536
  : RandomVariable (ExponentialVariableImpl (m, b))
692
{}
537
{}
693
double 
694
ExponentialVariable::GetSingleValue(double m, double b)
695
{
696
  return ExponentialVariableImpl::GetSingleValue (m, b);
697
}
698
538
699
//-----------------------------------------------------------------------------
539
//-----------------------------------------------------------------------------
700
//-----------------------------------------------------------------------------
540
//-----------------------------------------------------------------------------
 Lines 743-759    Link Here 
743
   */
583
   */
744
  virtual double GetValue();
584
  virtual double GetValue();
745
  virtual RandomVariableBase* Copy() const;
585
  virtual RandomVariableBase* Copy() const;
746
public:
586
747
  /**
748
   * \param m The mean value of the distribution from which the return value
749
   * is drawn.
750
   * \param s The shape parameter of the distribution from which the return
751
   * value is drawn.
752
   * \param b The upper bound to which to restrict return values
753
   * \return A random number from a Pareto distribution with mean m and shape
754
   * parameter s.
755
   */
756
  static double GetSingleValue(double m, double s, double b=0);
757
private:
587
private:
758
  double m_mean;  // Mean value of RV
588
  double m_mean;  // Mean value of RV
759
  double m_shape; // Shape parameter
589
  double m_shape; // Shape parameter
 Lines 778-797    Link Here 
778
608
779
double ParetoVariableImpl::GetValue()
609
double ParetoVariableImpl::GetValue()
780
{
610
{
781
  if(!RandomVariableBase::initialized)
782
  {
783
    RandomVariableBase::Initialize();
784
  }
785
  if(!m_generator)
611
  if(!m_generator)
786
  {
612
    {
787
    m_generator = new RngStream();
613
      m_generator = new RngStream();
788
    m_generator->InitializeStream();
614
    }
789
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
790
  }
791
  double scale = m_mean * ( m_shape - 1.0) / m_shape;
615
  double scale = m_mean * ( m_shape - 1.0) / m_shape;
792
  double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
616
  while(1)
793
  if (m_bound != 0 && r > m_bound) return m_bound;
617
    {
794
  return r;
618
      double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
619
      if (m_bound == 0 || r <= m_bound) return r;
620
      //otherwise, try again
621
    }
795
}
622
}
796
623
797
RandomVariableBase* ParetoVariableImpl::Copy() const
624
RandomVariableBase* ParetoVariableImpl::Copy() const
 Lines 799-818    Link Here 
799
  return new ParetoVariableImpl(*this);
626
  return new ParetoVariableImpl(*this);
800
}
627
}
801
628
802
double ParetoVariableImpl::GetSingleValue(double m, double s, double b/*=0*/)
803
{
804
  if(!RandomVariableBase::m_static_generator)
805
  {
806
    RandomVariableBase::Initialize(); // sets the static package seed
807
    RandomVariableBase::m_static_generator = new RngStream();
808
    RandomVariableBase::m_static_generator->InitializeStream();
809
  }
810
  double scale = m * ( s - 1.0) / s;
811
  double r = (scale * ( 1.0 / pow(m_static_generator->RandU01(), 1.0 / s)));
812
  if (b != 0 && r > b) return b;
813
  return r;
814
}
815
816
ParetoVariable::ParetoVariable ()
629
ParetoVariable::ParetoVariable ()
817
  : RandomVariable (ParetoVariableImpl ())
630
  : RandomVariable (ParetoVariableImpl ())
818
{}
631
{}
 Lines 825-835    Link Here 
825
ParetoVariable::ParetoVariable(double m, double s, double b)
638
ParetoVariable::ParetoVariable(double m, double s, double b)
826
  : RandomVariable (ParetoVariableImpl (m, s, b))
639
  : RandomVariable (ParetoVariableImpl (m, s, b))
827
{}
640
{}
828
double 
829
ParetoVariable::GetSingleValue(double m, double s, double b)
830
{
831
  return ParetoVariableImpl::GetSingleValue (m, s, b);
832
}
833
641
834
//-----------------------------------------------------------------------------
642
//-----------------------------------------------------------------------------
835
//-----------------------------------------------------------------------------
643
//-----------------------------------------------------------------------------
 Lines 879-892    Link Here 
879
   */
687
   */
880
  virtual double GetValue();
688
  virtual double GetValue();
881
  virtual RandomVariableBase* Copy(void) const;
689
  virtual RandomVariableBase* Copy(void) const;
882
public:
690
883
  /**
884
   * \param m Mean value for the distribution.
885
   * \param s Shape (alpha) parameter for the distribution.
886
   * \param b Upper limit on returned values
887
   * \return Random number from a distribution specified by m,s, and b
888
   */
889
  static double GetSingleValue(double m, double s, double b=0);
890
private:
691
private:
891
  double m_mean;  // Mean value of RV
692
  double m_mean;  // Mean value of RV
892
  double m_alpha; // Shape parameter
693
  double m_alpha; // Shape parameter
 Lines 906-925    Link Here 
906
707
907
double WeibullVariableImpl::GetValue()
708
double WeibullVariableImpl::GetValue()
908
{
709
{
909
  if(!RandomVariableBase::initialized)
910
  {
911
    RandomVariableBase::Initialize();
912
  }
913
  if(!m_generator)
710
  if(!m_generator)
914
  {
711
    {
915
    m_generator = new RngStream();
712
      m_generator = new RngStream();
916
    m_generator->InitializeStream();
713
    }
917
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
918
  }
919
  double exponent = 1.0 / m_alpha;
714
  double exponent = 1.0 / m_alpha;
920
  double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
715
  while(1)
921
  if (m_bound != 0 && r > m_bound) return m_bound;
716
    {
922
  return r;
717
      double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
718
      if (m_bound == 0 || r <= m_bound) return r;
719
      //otherwise, try again
720
    }
923
}
721
}
924
722
925
RandomVariableBase* WeibullVariableImpl::Copy() const
723
RandomVariableBase* WeibullVariableImpl::Copy() const
 Lines 927-946    Link Here 
927
  return new WeibullVariableImpl(*this);
725
  return new WeibullVariableImpl(*this);
928
}
726
}
929
727
930
double WeibullVariableImpl::GetSingleValue(double m, double s, double b/*=0*/)
931
{
932
  if(!RandomVariableBase::m_static_generator)
933
  {
934
    RandomVariableBase::Initialize(); // sets the static package seed
935
    RandomVariableBase::m_static_generator = new RngStream();
936
    RandomVariableBase::m_static_generator->InitializeStream();
937
  }
938
  double exponent = 1.0 / s;
939
  double r = m * pow( -log(m_static_generator->RandU01()), exponent);
940
  if (b != 0 && r > b) return b;
941
  return r;
942
}
943
944
WeibullVariable::WeibullVariable()
728
WeibullVariable::WeibullVariable()
945
  : RandomVariable (WeibullVariableImpl ())
729
  : RandomVariable (WeibullVariableImpl ())
946
{}
730
{}
 Lines 953-963    Link Here 
953
WeibullVariable::WeibullVariable(double m, double s, double b)
737
WeibullVariable::WeibullVariable(double m, double s, double b)
954
  : RandomVariable (WeibullVariableImpl (m, s, b))
738
  : RandomVariable (WeibullVariableImpl (m, s, b))
955
{}
739
{}
956
double 
740
957
WeibullVariable::GetSingleValue(double m, double s, double b)
958
{
959
  return WeibullVariableImpl::GetSingleValue (m, s, b);
960
}
961
//-----------------------------------------------------------------------------
741
//-----------------------------------------------------------------------------
962
//-----------------------------------------------------------------------------
742
//-----------------------------------------------------------------------------
963
// NormalVariableImpl methods
743
// NormalVariableImpl methods
 Lines 976-982    Link Here 
976
   * \brief Construct a normal random variable with specified mean and variance
756
   * \brief Construct a normal random variable with specified mean and variance
977
   * \param m Mean value
757
   * \param m Mean value
978
   * \param v Variance
758
   * \param v Variance
979
   * \param b Bound.  The NormalVariableImpl is bounded within +-bound.
759
   * \param b Bound.  The NormalVariableImpl is bounded within +-bound of the mean.
980
   */ 
760
   */ 
981
  NormalVariableImpl(double m, double v, double b = INFINITE_VALUE);
761
  NormalVariableImpl(double m, double v, double b = INFINITE_VALUE);
982
762
 Lines 987-1004    Link Here 
987
   */
767
   */
988
  virtual double GetValue();
768
  virtual double GetValue();
989
  virtual RandomVariableBase* Copy(void) const;
769
  virtual RandomVariableBase* Copy(void) const;
990
public:
770
991
  /**
771
  double GetMean (void) const;
992
   * \param m Mean value
772
  double GetVariance (void) const;
993
   * \param v Variance
773
  double GetBound (void) const;
994
   * \param b Bound.  The NormalVariableImpl is bounded within +-bound.
774
995
   * \return A random number from a distribution specified by m,v, and b.
996
   */ 
997
  static double GetSingleValue(double m, double v, double b = INFINITE_VALUE);
998
private:
775
private:
999
  double m_mean;      // Mean value of RV
776
  double m_mean;      // Mean value of RV
1000
  double m_variance;  // Mean value of RV
777
  double m_variance;  // Mean value of RV
1001
  double m_bound;     // Bound on value (absolute value)
778
  double m_bound;     // Bound on value's difference from the mean (absolute value)
1002
  bool   m_nextValid; // True if next valid
779
  bool   m_nextValid; // True if next valid
1003
  double m_next;      // The algorithm produces two values at a time
780
  double m_next;      // The algorithm produces two values at a time
1004
  static bool   m_static_nextValid;
781
  static bool   m_static_nextValid;
 Lines 1017-1036    Link Here 
1017
794
1018
NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c)
795
NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c)
1019
  : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance),
796
  : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance),
1020
    m_bound(c.m_bound), m_nextValid(false) { }
797
    m_bound(c.m_bound) { }
1021
798
1022
double NormalVariableImpl::GetValue()
799
double NormalVariableImpl::GetValue()
1023
{
800
{
1024
  if(!RandomVariableBase::initialized)
1025
  {
1026
    RandomVariableBase::Initialize();
1027
  }
1028
  if(!m_generator)
801
  if(!m_generator)
1029
  {
802
    {
1030
    m_generator = new RngStream();
803
      m_generator = new RngStream();
1031
    m_generator->InitializeStream();
804
    }
1032
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
1033
  }
1034
  if (m_nextValid)
805
  if (m_nextValid)
1035
    { // use previously generated
806
    { // use previously generated
1036
      m_nextValid = false;
807
      m_nextValid = false;
 Lines 1054-1068    Link Here 
1054
          double x1 = m_mean + v1 * y * sqrt(m_variance);
825
          double x1 = m_mean + v1 * y * sqrt(m_variance);
1055
          //if x1 is in bounds, return it
826
          //if x1 is in bounds, return it
1056
          if (fabs(x1-m_mean) <= m_bound)
827
          if (fabs(x1-m_mean) <= m_bound)
1057
          {
828
            {
1058
            return x1;
829
              return x1;
1059
          }
830
            }
1060
          //otherwise try and return m_next if it is valid
831
          //otherwise try and return m_next if it is valid
1061
          else if (m_nextValid)
832
          else if (m_nextValid)
1062
          {
833
	    {
1063
            m_nextValid = false;
834
	      m_nextValid = false;
1064
            return m_next;
835
	      return m_next;
1065
          }
836
	    }
1066
          //otherwise, just run this loop again
837
          //otherwise, just run this loop again
1067
        }
838
        }
1068
    }
839
    }
 Lines 1073-1121    Link Here 
1073
  return new NormalVariableImpl(*this);
844
  return new NormalVariableImpl(*this);
1074
}
845
}
1075
846
1076
double NormalVariableImpl::GetSingleValue(double m, double v, double b)
847
double
848
NormalVariableImpl::GetMean (void) const
1077
{
849
{
1078
  if(!RandomVariableBase::m_static_generator)
850
  return m_mean;
1079
  {
851
}
1080
    RandomVariableBase::Initialize(); // sets the static package seed
852
1081
    RandomVariableBase::m_static_generator = new RngStream();
853
double
1082
    RandomVariableBase::m_static_generator->InitializeStream();
854
NormalVariableImpl::GetVariance (void) const
1083
  }
855
{
1084
  if (m_static_nextValid)
856
  return m_variance;
1085
    { // use previously generated
857
}
1086
      m_static_nextValid = false;
858
1087
      return m_static_next;
859
double
1088
    }
860
NormalVariableImpl::GetBound (void) const
1089
  while(1)
861
{
1090
    { // See Simulation Modeling and Analysis p. 466 (Averill Law)
862
  return m_bound;
1091
      // for algorithm; basically a Box-Muller transform:
1092
      // http://en.wikipedia.org/wiki/Box-Muller_transform
1093
      double u1 = m_static_generator->RandU01();
1094
      double u2 = m_static_generator->RandU01();;
1095
      double v1 = 2 * u1 - 1;
1096
      double v2 = 2 * u2 - 1;
1097
      double w = v1 * v1 + v2 * v2;
1098
      if (w <= 1.0)
1099
        { // Got good pair
1100
          double y = sqrt((-2 * log(w))/w);
1101
          m_static_next = m + v2 * y * sqrt(v);
1102
          //if next is in bounds, it is valid
1103
          m_static_nextValid = fabs(m_static_next-m) <= b;;
1104
          double x1 = m + v1 * y * sqrt(v);
1105
          //if x1 is in bounds, return it
1106
          if (fabs(x1-m) <= b)
1107
          {
1108
            return x1;
1109
          }
1110
          //otherwise try and return m_next if it is valid
1111
          else if (m_static_nextValid)
1112
          {
1113
            m_static_nextValid = false;
1114
            return m_static_next;
1115
          }
1116
          //otherwise, just run this loop again
1117
        }
1118
    }
1119
}
863
}
1120
864
1121
NormalVariable::NormalVariable()
865
NormalVariable::NormalVariable()
 Lines 1127-1143    Link Here 
1127
NormalVariable::NormalVariable(double m, double v, double b)
871
NormalVariable::NormalVariable(double m, double v, double b)
1128
  : RandomVariable (NormalVariableImpl (m, v, b))
872
  : RandomVariable (NormalVariableImpl (m, v, b))
1129
{}
873
{}
1130
double 
1131
NormalVariable::GetSingleValue(double m, double v)
1132
{
1133
  return NormalVariableImpl::GetSingleValue (m, v);
1134
}
1135
double 
1136
NormalVariable::GetSingleValue(double m, double v, double b)
1137
{
1138
  return NormalVariableImpl::GetSingleValue (m, v, b);
1139
}
1140
1141
874
1142
//-----------------------------------------------------------------------------
875
//-----------------------------------------------------------------------------
1143
//-----------------------------------------------------------------------------
876
//-----------------------------------------------------------------------------
 Lines 1200-1220    Link Here 
1200
double EmpiricalVariableImpl::GetValue()
933
double EmpiricalVariableImpl::GetValue()
1201
{ // Return a value from the empirical distribution
934
{ // Return a value from the empirical distribution
1202
  // This code based (loosely) on code by Bruce Mah (Thanks Bruce!)
935
  // This code based (loosely) on code by Bruce Mah (Thanks Bruce!)
1203
  if(!RandomVariableBase::initialized)
1204
  {
1205
    RandomVariableBase::Initialize();
1206
  }
1207
  if(!m_generator)
936
  if(!m_generator)
1208
  {
937
    {
1209
    m_generator = new RngStream();
938
      m_generator = new RngStream();
1210
    m_generator->InitializeStream();
939
    }
1211
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
940
  if (emp.size() == 0) 
1212
  }
941
    {
1213
  if (emp.size() == 0) return 0.0; // HuH? No empirical data
942
      return 0.0; // HuH? No empirical data
1214
  if (!validated) Validate();      // Insure in non-decreasing
943
    }
944
  if (!validated) 
945
    {
946
      Validate();      // Insure in non-decreasing
947
    }
1215
  double r = m_generator->RandU01();
948
  double r = m_generator->RandU01();
1216
  if (r <= emp.front().cdf)return emp.front().value; // Less than first
949
  if (r <= emp.front().cdf)
1217
  if (r >= emp.back().cdf) return emp.back().value;  // Greater than last
950
    {
951
      return emp.front().value; // Less than first
952
    }
953
  if (r >= emp.back().cdf) 
954
    { 
955
      return emp.back().value;  // Greater than last
956
    }
1218
  // Binary search
957
  // Binary search
1219
  std::vector<ValueCDF>::size_type bottom = 0;
958
  std::vector<ValueCDF>::size_type bottom = 0;
1220
  std::vector<ValueCDF>::size_type top = emp.size() - 1;
959
  std::vector<ValueCDF>::size_type top = emp.size() - 1;
 Lines 1228-1235    Link Here 
1228
                             r);
967
                             r);
1229
        }
968
        }
1230
      // Not here, adjust bounds
969
      // Not here, adjust bounds
1231
      if (r < emp[c].cdf) top    = c - 1;
970
      if (r < emp[c].cdf)
1232
      else                bottom = c + 1;
971
        {
972
          top    = c - 1;
973
        }
974
      else
975
        {
976
          bottom = c + 1;
977
        }
1233
    }
978
    }
1234
}
979
}
1235
980
 Lines 1366-1372    Link Here 
1366
  
1111
  
1367
double DeterministicVariableImpl::GetValue()
1112
double DeterministicVariableImpl::GetValue()
1368
{
1113
{
1369
  if (next == count) next = 0;
1114
  if (next == count) 
1115
    {
1116
      next = 0;
1117
    }
1370
  return data[next++];
1118
  return data[next++];
1371
}
1119
}
1372
1120
 Lines 1395-1407    Link Here 
1395
   */
1143
   */
1396
  virtual double GetValue ();
1144
  virtual double GetValue ();
1397
  virtual RandomVariableBase* Copy(void) const;
1145
  virtual RandomVariableBase* Copy(void) const;
1398
public:
1146
1399
  /**
1400
   * \param mu mu parameter of the underlying normal distribution
1401
   * \param sigma sigma parameter of the underlying normal distribution
1402
   * \return A random number from the distribution specified by mu and sigma
1403
   */
1404
  static double GetSingleValue(double mu, double sigma);
1405
private:
1147
private:
1406
  double m_mu;
1148
  double m_mu;
1407
  double m_sigma;
1149
  double m_sigma;
 Lines 1447-1462    Link Here 
1447
double
1189
double
1448
LogNormalVariableImpl::GetValue ()
1190
LogNormalVariableImpl::GetValue ()
1449
{
1191
{
1450
  if(!RandomVariableBase::initialized)
1451
  {
1452
    RandomVariableBase::Initialize();
1453
  }
1454
  if(!m_generator)
1192
  if(!m_generator)
1455
  {
1193
    {
1456
    m_generator = new RngStream();
1194
      m_generator = new RngStream();
1457
    m_generator->InitializeStream();
1195
    }
1458
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
1459
  }
1460
  double u, v, r2, normal, z;
1196
  double u, v, r2, normal, z;
1461
1197
1462
  do
1198
  do
 Lines 1478-1519    Link Here 
1478
  return z;
1214
  return z;
1479
}
1215
}
1480
1216
1481
double LogNormalVariableImpl::GetSingleValue (double mu, double sigma)
1482
{
1483
  if(!RandomVariableBase::m_static_generator)
1484
  {
1485
    RandomVariableBase::Initialize(); // sets the static package seed
1486
    RandomVariableBase::m_static_generator = new RngStream();
1487
    RandomVariableBase::m_static_generator->InitializeStream();
1488
  }
1489
  double u, v, r2, normal, z;
1490
  do
1491
    {
1492
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
1493
      u = -1 + 2 * m_static_generator->RandU01 ();
1494
      v = -1 + 2 * m_static_generator->RandU01 ();
1495
1496
      /* see if it is in the unit circle */
1497
      r2 = u * u + v * v;
1498
    }
1499
  while (r2 > 1.0 || r2 == 0);
1500
1501
  normal = u * sqrt (-2.0 * log (r2) / r2);
1502
1503
  z =  exp (sigma * normal + mu);
1504
1505
  return z;
1506
}
1507
1508
LogNormalVariable::LogNormalVariable (double mu, double sigma)
1217
LogNormalVariable::LogNormalVariable (double mu, double sigma)
1509
  : RandomVariable (LogNormalVariableImpl (mu, sigma))
1218
  : RandomVariable (LogNormalVariableImpl (mu, sigma))
1510
{}
1219
{}
1511
double 
1512
LogNormalVariable::GetSingleValue(double mu, double sigma)
1513
{
1514
  return LogNormalVariableImpl::GetSingleValue (mu, sigma);
1515
}
1516
1517
1220
1518
//-----------------------------------------------------------------------------
1221
//-----------------------------------------------------------------------------
1519
//-----------------------------------------------------------------------------
1222
//-----------------------------------------------------------------------------
 Lines 1542-1555    Link Here 
1542
   */
1245
   */
1543
  virtual double GetValue();
1246
  virtual double GetValue();
1544
  virtual RandomVariableBase*  Copy(void) const;
1247
  virtual RandomVariableBase*  Copy(void) const;
1545
public:
1248
1546
  /**
1547
   * \param s Low end of the range
1548
   * \param l High end of the range
1549
   * \param mean mean of the distribution
1550
   * \return A triangularly distributed random number between s and l
1551
   */
1552
  static double GetSingleValue(double s, double l, double mean);
1553
private:
1249
private:
1554
  double m_min;
1250
  double m_min;
1555
  double m_max;
1251
  double m_max;
 Lines 1568-1588    Link Here 
1568
1264
1569
double TriangularVariableImpl::GetValue()
1265
double TriangularVariableImpl::GetValue()
1570
{
1266
{
1571
  if(!RandomVariableBase::initialized)
1572
  {
1573
    RandomVariableBase::Initialize();
1574
  }
1575
  if(!m_generator)
1267
  if(!m_generator)
1576
  {
1268
    {
1577
    m_generator = new RngStream();
1269
      m_generator = new RngStream();
1578
    m_generator->InitializeStream();
1270
    }
1579
    m_generator->ResetNthSubstream(RandomVariableBase::runNumber);
1580
  }
1581
  double u = m_generator->RandU01();
1271
  double u = m_generator->RandU01();
1582
  if(u <= (m_mode - m_min) / (m_max - m_min) )
1272
  if(u <= (m_mode - m_min) / (m_max - m_min) )
1583
    return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) );
1273
    {
1274
      return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) );
1275
    }
1584
  else
1276
  else
1585
    return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) );
1277
    {
1278
      return m_max - sqrt( (1-u) * (m_max - m_min) * (m_max - m_mode) );
1279
    }
1586
}
1280
}
1587
1281
1588
RandomVariableBase* TriangularVariableImpl::Copy() const
1282
RandomVariableBase* TriangularVariableImpl::Copy() const
 Lines 1590-1623    Link Here 
1590
  return new TriangularVariableImpl(*this);
1284
  return new TriangularVariableImpl(*this);
1591
}
1285
}
1592
1286
1593
double TriangularVariableImpl::GetSingleValue(double s, double l, double mean)
1594
{
1595
  if(!RandomVariableBase::m_static_generator)
1596
  {
1597
    RandomVariableBase::Initialize(); // sets the static package seed
1598
    RandomVariableBase::m_static_generator = new RngStream();
1599
    RandomVariableBase::m_static_generator->InitializeStream();
1600
  }
1601
  double mode = 3.0*mean-s-l;
1602
  double u = m_static_generator->RandU01();
1603
  if(u <= (mode - s) / (l - s) )
1604
    return s + sqrt(u * (l - s) * (mode - s) );
1605
  else
1606
    return l - sqrt( (1-u) * (l - s) * (l - mode) );
1607
}
1608
1609
TriangularVariable::TriangularVariable()
1287
TriangularVariable::TriangularVariable()
1610
  : RandomVariable (TriangularVariableImpl ())
1288
  : RandomVariable (TriangularVariableImpl ())
1611
{}
1289
{}
1612
TriangularVariable::TriangularVariable(double s, double l, double mean)
1290
TriangularVariable::TriangularVariable(double s, double l, double mean)
1613
  : RandomVariable (TriangularVariableImpl (s,l,mean))
1291
  : RandomVariable (TriangularVariableImpl (s,l,mean))
1614
{}
1292
{}
1615
double 
1293
1616
TriangularVariable::GetSingleValue(double s, double l, double mean)
1294
//-----------------------------------------------------------------------------
1295
//-----------------------------------------------------------------------------
1296
// ZipfVariableImpl
1297
class ZipfVariableImpl : public RandomVariableBase { 
1298
public:
1299
  /**
1300
   * \param n the number of possible items
1301
   * \param alpha the alpha parameter
1302
   */
1303
  ZipfVariableImpl (long n, double alpha);
1304
1305
  /**
1306
   * \A zipf variable with N=1 and alpha=0
1307
   */
1308
  ZipfVariableImpl ();
1309
1310
  /**
1311
   * \return A random value from this distribution
1312
   */
1313
  virtual double GetValue ();
1314
  virtual RandomVariableBase* Copy(void) const;
1315
1316
private:
1317
  long m_n;
1318
  double m_alpha;
1319
  double m_c; //the normalization constant
1320
};
1321
1322
1323
RandomVariableBase* ZipfVariableImpl::Copy () const
1617
{
1324
{
1618
  return TriangularVariableImpl::GetSingleValue (s,l,mean);
1325
  return new ZipfVariableImpl (m_n, m_alpha);
1619
}
1326
}
1620
1327
1328
ZipfVariableImpl::ZipfVariableImpl ()
1329
    :m_n(1), m_alpha(0), m_c(1)
1330
{
1331
}
1332
1333
1334
ZipfVariableImpl::ZipfVariableImpl (long n, double alpha)
1335
    :m_n(n), m_alpha(alpha), m_c(0)
1336
{
1337
  //calculate the normalization constant c
1338
  for(int i=1;i<=n;i++)
1339
    m_c+=(1.0/pow((double)i,alpha));
1340
  m_c=1.0/m_c;
1341
}
1342
1343
double
1344
ZipfVariableImpl::GetValue ()
1345
{
1346
  if(!m_generator)
1347
  {
1348
    m_generator = new RngStream();
1349
  }
1350
1351
  double u = m_generator->RandU01();
1352
  double sum_prob=0,zipf_value=0;
1353
  for(int i=1;i<=m_n;i++)
1354
  {
1355
    sum_prob+=m_c/pow((double)i,m_alpha);
1356
    if(sum_prob>u)
1357
    {
1358
      zipf_value=i;
1359
      break;
1360
    }
1361
  }
1362
  return zipf_value;
1363
}
1364
1365
ZipfVariable::ZipfVariable ()
1366
  : RandomVariable (ZipfVariableImpl ())
1367
{}
1368
1369
ZipfVariable::ZipfVariable (long n, double alpha)
1370
  : RandomVariable (ZipfVariableImpl (n, alpha))
1371
{}
1372
1621
1373
1622
std::ostream &operator << (std::ostream &os, const RandomVariable &var)
1374
std::ostream &operator << (std::ostream &os, const RandomVariable &var)
1623
{
1375
{
 Lines 1634-1639    Link Here 
1634
      os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax ();
1386
      os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax ();
1635
      return os;
1387
      return os;
1636
    }
1388
    }
1389
  NormalVariableImpl *normal = dynamic_cast<NormalVariableImpl *> (base);
1390
  if (normal != 0)
1391
    {
1392
      os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance ();
1393
      double bound = normal->GetBound ();
1394
      if (bound != NormalVariableImpl::INFINITE_VALUE)
1395
        {
1396
          os << ":" << bound;
1397
        }
1398
      return os;
1399
    }
1637
  // XXX: support other distributions
1400
  // XXX: support other distributions
1638
  os.setstate (std::ios_base::badbit);
1401
  os.setstate (std::ios_base::badbit);
1639
  return os;
1402
  return os;
 Lines 1679-1684    Link Here 
1679
          var = UniformVariable (a, b);
1442
          var = UniformVariable (a, b);
1680
        }
1443
        }
1681
    }
1444
    }
1445
  else if (type == "Normal")
1446
    {
1447
      if (value.size () == 0)
1448
        {
1449
          var = NormalVariable ();
1450
        }
1451
      else
1452
        {
1453
          tmp = value.find (":");
1454
          if (tmp == value.npos)
1455
            {
1456
              NS_FATAL_ERROR ("bad Normal value: " << value);
1457
            }
1458
          std::string::size_type tmp2;
1459
          std::string sub = value.substr (tmp + 1, value.npos);
1460
          tmp2 = sub.find (":");
1461
          if (tmp2 == value.npos)
1462
            {
1463
              istringstream issA (value.substr (0, tmp));
1464
              istringstream issB (sub);
1465
              double a, b;
1466
              issA >> a;
1467
              issB >> b;
1468
              var = NormalVariable (a, b);
1469
            }
1470
          else
1471
            {
1472
              istringstream issA (value.substr (0, tmp));
1473
              istringstream issB (sub.substr (0, tmp2));
1474
              istringstream issC (sub.substr (tmp2 + 1, value.npos));
1475
              double a, b, c;
1476
              issA >> a;
1477
              issB >> b;
1478
              issC >> c;
1479
              var = NormalVariable (a, b, c);
1480
            }
1481
        }
1482
    }
1682
  else
1483
  else
1683
    {
1484
    {
1684
      NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << type);
1485
      NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << type);
 Lines 1691-1696    Link Here 
1691
1492
1692
}//namespace ns3
1493
}//namespace ns3
1693
1494
1495
 
1694
1496
1695
#ifdef RUN_SELF_TESTS
1497
#ifdef RUN_SELF_TESTS
1696
#include "test.h"
1498
#include "test.h"
 Lines 1747-1793    Link Here 
1747
        }
1549
        }
1748
    }
1550
    }
1749
1551
1750
    // Test GetSingleValue
1751
    {
1752
      vector<double> samples;
1753
      const int NSAMPLES = 10000;
1754
      double sum = 0;
1755
      for (int n = NSAMPLES; n; --n)
1756
        {
1757
          double value = LogNormalVariable::GetSingleValue (mu, sigma);
1758
          sum += value;
1759
          samples.push_back (value);
1760
        }
1761
      double obtained_mean = sum / NSAMPLES;
1762
      sum = 0;
1763
      for (vector<double>::iterator iter = samples.begin (); iter != samples.end (); iter++)
1764
        {
1765
          double tmp = (*iter - obtained_mean);
1766
          sum += tmp*tmp;
1767
        }
1768
      double obtained_stddev = sqrt (sum / (NSAMPLES - 1));
1769
1770
      if (not (obtained_mean/desired_mean > 0.90 and obtained_mean/desired_mean < 1.10))
1771
        {
1772
          result = false;
1773
          Failure () << "Obtained LogNormalVariable::GetSingleValue mean value " << obtained_mean
1774
                     << ", expected " << desired_mean << std::endl;
1775
        }
1776
1777
      if (not (obtained_stddev/desired_stddev > 0.90 and obtained_stddev/desired_stddev < 1.10))
1778
        {
1779
          result = false;
1780
          Failure () << "Obtained LogNormalVariable::GetSingleValue stddev value " << obtained_stddev <<
1781
            ", expected " << desired_stddev << std::endl;
1782
        }
1783
    }
1784
1785
    // Test attribute serialization
1552
    // Test attribute serialization
1786
    {
1553
    {
1787
      RandomVariableValue val;
1554
      {
1788
      val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ());
1555
        RandomVariableValue val;
1789
      RandomVariable rng = val.Get ();
1556
        val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ());
1790
      NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2");
1557
        RandomVariable rng = val.Get ();
1558
        NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Uniform:0.1:0.2");
1559
      }
1560
      {
1561
        RandomVariableValue val;
1562
        val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ());
1563
        RandomVariable rng = val.Get ();
1564
        NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2");
1565
      }
1566
      {
1567
        RandomVariableValue val;
1568
        val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker ());
1569
        RandomVariable rng = val.Get ();
1570
        NS_TEST_ASSERT_EQUAL (val.SerializeToString (MakeRandomVariableChecker ()), "Normal:0.1:0.2:0.15");
1571
      }
1791
    }
1572
    }
1792
1573
1793
    return result;
1574
    return result;
(-)a/src/core/random-variable.h (-144 / +108 lines)
 Lines 16-25    Link Here 
16
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
16
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17
//
17
//
18
// Author: Rajib Bhattacharjea<raj.b@gatech.edu>
18
// Author: Rajib Bhattacharjea<raj.b@gatech.edu>
19
// Author: Hadi Arbabi<marbabi@cs.odu.edu>
19
//
20
//
20
21
21
#ifndef __random_variable_h__
22
#ifndef __random_variable_h
22
#define __random_variable_h__
23
#define __random_variable_h
23
24
24
#include <vector>
25
#include <vector>
25
#include <algorithm>
26
#include <algorithm>
 Lines 39-44    Link Here 
39
40
40
class RandomVariableBase;
41
class RandomVariableBase;
41
42
43
class SeedManager
44
{
45
public:
46
	  
47
  /**
48
   * \brief set the seed
49
   * it will duplicate the seed value 6 times
50
   * \code
51
   * SeedManger::SetSeed(15);
52
   * UniformVariable x(2,3);     //these will give the same output everytime
53
   * ExponentialVariable y(120); //as long as the seed stays the same
54
   * \endcode
55
   * \param seed
56
   *
57
   * Note, while the underlying RNG takes six integer values as a seed;
58
   * it is sufficient to set these all to the same integer, so we provide
59
   * a simpler interface here that just takes one integer.
60
   */ 
61
  static void SetSeed (uint32_t seed);
62
 
63
  /**
64
   * \brief Get the seed value
65
   * \return the seed value
66
   *
67
   * Note:  returns the first of the six seed values used in the underlying RNG
68
   */
69
   static uint32_t GetSeed ();
70
 
71
   /**
72
    * \brief Set the run number of simulation
73
    *
74
    * \code
75
    * SeedManager::SetSeed(12);
76
    * int N = atol(argv[1]); //read in run number from command line
77
    * SeedManager::SetRun(N);
78
    * UniformVariable x(0,10);
79
    * ExponentialVariable y(2902);
80
    * \endcode
81
    * In this example, N could successivly be equal to 1,2,3, etc. and the user
82
    * would continue to get independent runs out of the single simulation.  For
83
    * this simple example, the following might work:
84
    * \code
85
    * ./simulation 0
86
    * ...Results for run 0:...
87
    *
88
    * ./simulation 1
89
    * ...Results for run 1:...
90
    * \endcode
91
    */
92
  static void SetRun (uint32_t run);
93
  /**
94
   * \returns the current run number
95
   * @sa SetRun
96
   */
97
  static uint32_t GetRun (void);
98
  
99
  /**
100
   * \brief Check if seed value is valid if wanted to be used as seed
101
   * \return true if valid and false if invalid
102
   */
103
  static bool CheckSeed (uint32_t seed);
104
};
105
106
42
/**
107
/**
43
 * \brief The basic RNG for NS-3.
108
 * \brief The basic RNG for NS-3.
44
 * \ingroup randomvariable
109
 * \ingroup randomvariable
 Lines 73-167    Link Here 
73
   * \return  Integer cast of ::GetValue()
138
   * \return  Integer cast of ::GetValue()
74
   */
139
   */
75
  uint32_t GetInteger (void) const;
140
  uint32_t GetInteger (void) const;
76
  
77
  /**
78
   * \brief Get the internal state of the RNG
79
   *
80
   * This function is for power users who understand the inner workings
81
   * of the underlying RngStream method used.  It returns the internal
82
   * state of the RNG via the input parameter.
83
   * \param seed Output parameter; gets overwritten with the internal state of
84
   * of the RNG.
85
   */
86
  void GetSeed(uint32_t seed[6]) const;
87
  
88
  /**
89
   * \brief Set seeding behavior
90
   * 
91
   * Specify whether the POSIX device /dev/random is to
92
   * be used for seeding.  When this is used, the underlying
93
   * generator is seeded with data from /dev/random instead of
94
   * being seeded based upon the time of day.  For this to be effective,
95
   * it must be called before the creation of the first instance of a 
96
   * RandomVariable or subclass.  Example:
97
   * \code
98
   * RandomVariable::UseDevRandom();
99
   * UniformVariable x(2,3);  //these are seeded randomly
100
   * ExponentialVariable y(120); //etc
101
   * \endcode
102
   * \param udr True if /dev/random desired.
103
   */
104
  static  void UseDevRandom(bool udr = true);
105
106
   /**
107
   * \brief Use the global seed to force precisely reproducible results.
108
   *
109
   * It is often desirable to create a simulation that uses random
110
   * numbers, while at the same time is completely reproducible.
111
   * Specifying this set of six random seeds initializes the
112
   * random number generator with the specified seed.
113
   * Once this is set, all generators will produce fixed output
114
   * from run to run.  This is because each time a new generator is created,
115
   * the underlying RngStream deterministically creates a new seed based upon
116
   * the old one, hence a "stream" of RNGs.  Example:
117
   * \code
118
   * RandomVariable::UseGlobalSeed(...);
119
   * UniformVariable x(2,3);     //these will give the same output everytime
120
   * ExponentialVariable y(120); //as long as the seed stays the same
121
   * \endcode
122
   * \param s0
123
   * \param s1
124
   * \param s2
125
   * \param s3
126
   * \param s4
127
   * \param s5
128
   * \return True if seed is valid.
129
   */ 
130
  static void UseGlobalSeed(uint32_t s0, uint32_t s1, uint32_t s2, 
131
                            uint32_t s3, uint32_t s4, uint32_t s5);
132
  
133
  /**
134
   * \brief Set the run number of this simulation
135
   *
136
   * These RNGs have the ability to give independent sets of trials for a fixed
137
   * global seed.  For example, suppose one sets up a simulation with
138
   * RandomVariables with a given global seed.  Suppose the user wanted to
139
   * retry the same simulation with different random values for validity,
140
   * statistical rigor, etc.  The user could either change the global seed and
141
   * re-run the simulation, or could use this facility to increment all of the
142
   * RNGs to a next substream state.  This predictably advances the internal
143
   * state of all RandomVariables n steps.  This should be called immediately
144
   * after the global seed is set, and before the creation of any
145
   * RandomVariables.  For example:
146
   * \code
147
   * RandomVariable::UseGlobalSeed(1,2,3,4,5,6);
148
   * int N = atol(argv[1]); //read in run number from command line
149
   * RandomVariable::SetRunNumber(N);
150
   * UniformVariable x(0,10);
151
   * ExponentialVariable y(2902);
152
   * \endcode
153
   * In this example, N could successivly be equal to 1,2,3, etc. and the user
154
   * would continue to get independent runs out of the single simulation.  For
155
   * this simple example, the following might work:
156
   * \code
157
   * ./simulation 0
158
   * ...Results for run 0:...
159
   *
160
   * ./simulation 1
161
   * ...Results for run 1:...
162
   * \endcode
163
   */
164
  static void SetRunNumber(uint32_t n);
165
141
166
private:
142
private:
167
  friend std::ostream &operator << (std::ostream &os, const RandomVariable &var);
143
  friend std::ostream &operator << (std::ostream &os, const RandomVariable &var);
 Lines 204-216    Link Here 
204
   * \param l High end of the range
180
   * \param l High end of the range
205
   */
181
   */
206
  UniformVariable(double s, double l);
182
  UniformVariable(double s, double l);
207
public:
183
208
  /**
184
  /**
185
  * \brief call RandomVariable::GetValue
186
  * \return A floating point random value
187
  *
188
  * Note: we have to re-implement this method here because the method is 
189
  * overloaded below for the two-argument variant and the c++ name resolution
190
  * rules don't work well with overloads split between parent and child 
191
  * classes.
192
  */
193
  double GetValue (void) const;
194
  
195
  /**
196
  * \brief Returns a random double with the specified range
197
  * \param s Low end of the range
198
  * \param l High end of the range
199
  * \return A floating point random value
200
  */
201
  double GetValue(double s, double l);
202
203
  /**
204
   * \brief Returns a random unsigned integer from the interval [s,l] including both ends.
209
   * \param s Low end of the range
205
   * \param s Low end of the range
210
   * \param l High end of the range
206
   * \param l High end of the range
211
   * \return A uniformly distributed random number between s and l
207
   * \return A random unsigned integer value.
212
   */
208
   */
213
  static double GetSingleValue(double s, double l);
209
  uint32_t GetInteger (uint32_t s, uint32_t l);
214
};
210
};
215
211
216
/**
212
/**
 Lines 334-345    Link Here 
334
   */
330
   */
335
  ExponentialVariable(double m, double b);
331
  ExponentialVariable(double m, double b);
336
332
337
  /**
338
   * \param m The mean of the distribution from which the return value is drawn
339
   * \param b The upper bound value desired, beyond which values get clipped
340
   * \return A random number from an exponential distribution with mean m
341
   */
342
  static double GetSingleValue(double m, double b=0);
343
};
333
};
344
334
345
/**
335
/**
 Lines 402-417    Link Here 
402
   */
392
   */
403
  ParetoVariable(double m, double s, double b);
393
  ParetoVariable(double m, double s, double b);
404
394
405
  /**
406
   * \param m The mean value of the distribution from which the return value
407
   * is drawn.
408
   * \param s The shape parameter of the distribution from which the return
409
   * value is drawn.
410
   * \param b The upper bound to which to restrict return values
411
   * \return A random number from a Pareto distribution with mean m and shape
412
   * parameter s.
413
   */
414
  static double GetSingleValue(double m, double s, double b=0);
415
};
395
};
416
396
417
/**
397
/**
 Lines 464-476    Link Here 
464
   * \param b Upper limit on returned values
444
   * \param b Upper limit on returned values
465
   */
445
   */
466
  WeibullVariable(double m, double s, double b);
446
  WeibullVariable(double m, double s, double b);
467
  /**
447
468
   * \param m Mean value for the distribution.
469
   * \param s Shape (alpha) parameter for the distribution.
470
   * \param b Upper limit on returned values
471
   * \return Random number from a distribution specified by m,s, and b
472
   */
473
  static double GetSingleValue(double m, double s, double b=0);
474
};
448
};
475
449
476
/**
450
/**
 Lines 511-532    Link Here 
511
   * [mean-bound,mean+bound]
485
   * [mean-bound,mean+bound]
512
   */ 
486
   */ 
513
  NormalVariable(double m, double v, double b);
487
  NormalVariable(double m, double v, double b);
514
515
  /**
516
   * \param m Mean value
517
   * \param v Variance
518
   * \return A random number from a distribution specified by m, and v.
519
   */ 
520
  static double GetSingleValue(double m, double v);
521
522
  /**
523
   * \param m Mean value
524
   * \param v Variance
525
   * \param b Bound.  The NormalVariable is bounded symetrically about the mean
526
   * [mean-bound,mean+bound]
527
   * \return A random number from a distribution specified by m,v, and b.
528
   */ 
529
  static double GetSingleValue(double m, double v, double b);
530
};
488
};
531
489
532
/**
490
/**
 Lines 633-645    Link Here 
633
   * \param sigma sigma parameter of the lognormal distribution
591
   * \param sigma sigma parameter of the lognormal distribution
634
   */
592
   */
635
  LogNormalVariable (double mu, double sigma);
593
  LogNormalVariable (double mu, double sigma);
594
};
636
595
596
/**
597
 * \brief Zipf Distributed random var (between 1 and n included)
598
 * \ingroup randomvariable
599
 *
600
 */
601
class ZipfVariable : public RandomVariable 
602
{
603
public:
637
  /**
604
  /**
638
   * \param mu mu parameter of the underlying normal distribution
605
   * \param n the number of possible items
639
   * \param sigma sigma parameter of the underlying normal distribution
606
   * \param alpha the alpha parameter
640
   * \return A random number from the distribution specified by mu and sigma
641
   */
607
   */
642
  static double GetSingleValue(double mu, double sigma);
608
  ZipfVariable (long n, double alpha);
609
  /**
610
   * A zipf variable with N=1 and alpha=0
611
   */
612
  ZipfVariable ();
643
};
613
};
644
614
645
/**
615
/**
 Lines 666-678    Link Here 
666
   * \param mean mean of the distribution
636
   * \param mean mean of the distribution
667
   */
637
   */
668
  TriangularVariable(double s, double l, double mean);
638
  TriangularVariable(double s, double l, double mean);
669
  /**
639
670
   * \param s Low end of the range
671
   * \param l High end of the range
672
   * \param mean mean of the distribution
673
   * \return A triangularly distributed random number between s and l
674
   */
675
  static double GetSingleValue(double s, double l, double mean);
676
};
640
};
677
641
678
std::ostream &operator << (std::ostream &os, const RandomVariable &var);
642
std::ostream &operator << (std::ostream &os, const RandomVariable &var);

Return to bug 578