|
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> |
|
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" |
|
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 |
|
|
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() |
|
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) |
|
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 |
|
|
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; |
|
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 |
//----------------------------------------------------------------------------- |
|
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) |
|
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 ()) |
|
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 |
//----------------------------------------------------------------------------- |
|
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 |
|
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 |
|
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 |
{} |
|
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 |
//----------------------------------------------------------------------------- |
|
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 |
|
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 |
|
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 |
{} |
|
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 |
|
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 |
|
|
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; |
|
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; |
|
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 |
} |
|
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() |
|
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 |
//----------------------------------------------------------------------------- |
|
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; |
|
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 |
|
|
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 |
|
|
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; |
|
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 |
|
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 |
//----------------------------------------------------------------------------- |
|
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; |
|
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 |
|
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 |
{ |
|
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; |
|
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); |
|
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" |
|
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; |