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