Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(1014)

Side by Side Diff: src/core/random-variable.cc

Issue 193060: random-variable - doc updated and new Zeta distribution (Closed)
Patch Set: fixed last stuff as per Michele Weigle comments Created 14 years, 12 months ago
Left:
Right:
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View unified diff | Download patch
« no previous file with comments | « src/core/random-variable.h ('k') | no next file » | no next file with comments »
Toggle Intra-line Diffs ('i') | Expand Comments ('e') | Collapse Comments ('c') | Show Comments Hide Comments ('s')
OLDNEW
1 /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */ 1 /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
2 // 2 //
3 // Copyright (c) 2006 Georgia Tech Research Corporation 3 // Copyright (c) 2006 Georgia Tech Research Corporation
4 // 4 //
5 // This program is free software; you can redistribute it and/or modify 5 // This program is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License version 2 as 6 // it under the terms of the GNU General Public License version 2 as
7 // published by the Free Software Foundation; 7 // published by the Free Software Foundation;
8 // 8 //
9 // This program is distributed in the hope that it will be useful, 9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of 10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details. 12 // GNU General Public License for more details.
13 // 13 //
14 // You should have received a copy of the GNU General Public License 14 // You should have received a copy of the GNU General Public License
15 // along with this program; if not, write to the Free Software 15 // along with this program; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 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 // Author: Hadi Arbabi<marbabi@cs.odu.edu>
20 // 20 //
21 21
22 #include <iostream> 22 #include <iostream>
23 23
24 #include <math.h> 24 #include <math.h>
25 #include <stdlib.h> 25 #include <stdlib.h>
26 #include <sys/time.h>» » » // for gettimeofday 26 #include <sys/time.h> // for gettimeofday
27 #include <unistd.h> 27 #include <unistd.h>
28 #include <iostream> 28 #include <iostream>
29 #include <sys/types.h> 29 #include <sys/types.h>
30 #include <sys/stat.h> 30 #include <sys/stat.h>
31 #include <fcntl.h> 31 #include <fcntl.h>
32 #include <sstream> 32 #include <sstream>
33 #include <vector> 33 #include <vector>
34 34
35 #include "test.h" 35 #include "test.h"
36 #include "assert.h" 36 #include "assert.h"
37 #include "config.h" 37 #include "config.h"
38 #include "integer.h" 38 #include "integer.h"
39 #include "random-variable.h" 39 #include "random-variable.h"
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 {
144 return *this; 147 return *this;
145 } 148 }
146 delete m_variable; 149 delete m_variable;
147 m_variable = o.m_variable->Copy (); 150 m_variable = o.m_variable->Copy ();
148 return *this; 151 return *this;
149 } 152 }
150 RandomVariable::~RandomVariable() 153 RandomVariable::~RandomVariable ()
151 { 154 {
152 delete m_variable; 155 delete m_variable;
153 } 156 }
154 double 157 double
155 RandomVariable::GetValue (void) const 158 RandomVariable::GetValue (void) const
156 { 159 {
157 return m_variable->GetValue (); 160 return m_variable->GetValue ();
158 } 161 }
159 162
160 uint32_t 163 uint32_t
161 RandomVariable::GetInteger (void) const 164 RandomVariable::GetInteger (void) const
162 { 165 {
163 return m_variable->GetInteger (); 166 return m_variable->GetInteger ();
164 } 167 }
165 168
166 RandomVariableBase * 169 RandomVariableBase *
167 RandomVariable::Peek (void) const 170 RandomVariable::Peek (void) const
168 { 171 {
169 return m_variable; 172 return m_variable;
170 } 173 }
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 *
372 * The four parameters define the sequence. For example 401 * The four parameters define the sequence. For example
373 * SequentialVariableImpl(0,5,1,2) creates a RNG that has the sequence 402 * SequentialVariableImpl(0,5,1,2) creates a RNG that has the sequence
374 * 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0 ... 403 * 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 0, 0 ...
375 * \param f First value of the sequence. 404 * \param f First value of the sequence.
376 * \param l One more than the last value of the sequence. 405 * \param l One more than the last value of the sequence.
377 * \param i Increment between sequence values 406 * \param i Increment between sequence values
378 * \param c Number of times each member of the sequence is repeated 407 * \param c Number of times each member of the sequence is repeated
379 */ 408 */
380 SequentialVariableImpl(double f, double l, double i = 1, uint32_t c = 1); 409 SequentialVariableImpl (double f, double l, double i = 1, uint32_t c = 1);
381 410
382 /** 411 /**
383 * \brief Constructor for the SequentialVariableImpl RNG. 412 * \brief Constructor for the SequentialVariableImpl RNG.
384 * 413 *
385 * Differs from the first only in that the increment parameter is a 414 * Differs from the first only in that the increment parameter is a
386 * random variable 415 * random variable
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, uin t32_t c) 440 SequentialVariableImpl::SequentialVariableImpl (double f, double l, double i, ui nt32_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)),
444 m_consecutive (c),
445 m_current (f),
446 m_currentConsecutive (0)
447 {
448 }
415 449
416 SequentialVariableImpl::SequentialVariableImpl(double f, double l, const RandomV ariable& i, uint32_t c) 450 SequentialVariableImpl::SequentialVariableImpl (double f, double l, const Random Variable& i, uint32_t c)
417 : m_min(f), m_max(l), m_increment(i), m_consecutive(c), 451 : m_min (f),
418 m_current(f), m_currentConsecutive(0) 452 m_max (l),
419 {} 453 m_increment (i),
454 m_consecutive (c),
455 m_current (f),
456 m_currentConsecutive (0)
457 {
458 }
420 459
421 SequentialVariableImpl::SequentialVariableImpl(const SequentialVariableImpl& c) 460 SequentialVariableImpl::SequentialVariableImpl (const SequentialVariableImpl& c)
422 : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max), 461 : RandomVariableBase (c),
423 m_increment(c.m_increment), m_consecutive(c.m_consecutive), 462 m_min (c.m_min),
424 m_current(c.m_current), m_currentConsecutive(c.m_currentConsecutive) 463 m_max (c.m_max),
425 {} 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 }
426 470
427 SequentialVariableImpl::~SequentialVariableImpl() 471 SequentialVariableImpl::~SequentialVariableImpl ()
428 {} 472 {
473 }
429 474
430 double SequentialVariableImpl::GetValue() 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 {
502 }
454 503
455 //----------------------------------------------------------------------------- 504 // -----------------------------------------------------------------------------
456 //----------------------------------------------------------------------------- 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 ··
505 ExponentialVariableImpl::ExponentialVariableImpl(double m, double b)·
506 : m_mean(m), m_bound(b) { }
507 ··
508 ExponentialVariableImpl::ExponentialVariableImpl(const ExponentialVariableImpl& c)·
509 : RandomVariableBase(c), m_mean(c.m_mean), m_bound(c.m_bound) { }
510 554
511 double ExponentialVariableImpl::GetValue() 555 ExponentialVariableImpl::ExponentialVariableImpl (double m)
556 : m_mean (m),
557 m_bound (0)
512 { 558 {
513 if(!m_generator) 559 }
560
561 ExponentialVariableImpl::ExponentialVariableImpl (double m, double b)
562 : m_mean (m),
563 m_bound (b)
564 {
565 }
566
567 ExponentialVariableImpl::ExponentialVariableImpl (const ExponentialVariableImpl& c)
568 : RandomVariableBase (c),
569 m_mean (c.m_mean),
570 m_bound (c.m_bound)
571 {
572 }
573
574 double ExponentialVariableImpl::GetValue ()
575 {
576 if (!m_generator)
514 { 577 {
515 m_generator = new RngStream(); 578 m_generator = new RngStream ();
516 } 579 }
517 while(1) 580 while (1)
518 { 581 {
519 double r = -m_mean*log(m_generator->RandU01()); 582 double r = -m_mean*log (m_generator->RandU01 ());
520 if (m_bound == 0 || r <= m_bound) return r; 583 if (m_bound == 0 || r <= m_bound)
521 //otherwise, try again 584 {
585 return r;
586 }
587 // otherwise, try again
522 } 588 }
523 } 589 }
524 590
525 RandomVariableBase* ExponentialVariableImpl::Copy() const 591 RandomVariableBase* ExponentialVariableImpl::Copy () const
526 { 592 {
527 return new ExponentialVariableImpl(*this); 593 return new ExponentialVariableImpl (*this);
528 } 594 }
529 595
530 ExponentialVariable::ExponentialVariable() 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 {
607 }
539 608
540 //----------------------------------------------------------------------------- 609 // -----------------------------------------------------------------------------
541 //----------------------------------------------------------------------------- 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 *
570 * Since pareto distributions can theoretically return unbounded values, 640 * Since pareto distributions can theoretically return unbounded values,
571 * it is sometimes useful to specify a fixed upper limit. Note however 641 * it is sometimes useful to specify a fixed upper limit. Note however
572 * when the upper limit is specified, the true mean of the distribution 642 * when the upper limit is specified, the true mean of the distribution
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),
666 m_shape (1.5),
667 m_bound (0)
668 {
669 }
596 670
597 ParetoVariableImpl::ParetoVariableImpl(double m) 671 ParetoVariableImpl::ParetoVariableImpl (double m)
598 : m_mean(m), m_shape(1.5), m_bound(0) { } 672 : m_mean (m),
673 m_shape (1.5),
674 m_bound (0)
675 {
676 }
599 677
600 ParetoVariableImpl::ParetoVariableImpl(double m, double s) 678 ParetoVariableImpl::ParetoVariableImpl (double m, double s)
601 : m_mean(m), m_shape(s), m_bound(0) { } 679 : m_mean (m),
680 m_shape (s),
681 m_bound (0)
682 {
683 }
602 684
603 ParetoVariableImpl::ParetoVariableImpl(double m, double s, double b) 685 ParetoVariableImpl::ParetoVariableImpl (double m, double s, double b)
604 : m_mean(m), m_shape(s), m_bound(b) { } 686 : m_mean (m),
687 m_shape (s),
688 m_bound (b)
689 {
690 }
605 691
606 ParetoVariableImpl::ParetoVariableImpl(const ParetoVariableImpl& c) 692 ParetoVariableImpl::ParetoVariableImpl (const ParetoVariableImpl& c)
607 : RandomVariableBase(c), m_mean(c.m_mean), m_shape(c.m_shape), 693 : RandomVariableBase (c),
608 m_bound(c.m_bound) { } 694 m_mean (c.m_mean),
695 m_shape (c.m_shape),
696 m_bound (c.m_bound)
697 {
698 }
609 699
610 double ParetoVariableImpl::GetValue() 700 double ParetoVariableImpl::GetValue ()
611 { 701 {
612 if(!m_generator) 702 if (!m_generator)
613 { 703 {
614 m_generator = new RngStream(); 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 {
712 return r;
713 }
714 // otherwise, try again
622 } 715 }
623 } 716 }
624 717
625 RandomVariableBase* ParetoVariableImpl::Copy() const 718 RandomVariableBase* ParetoVariableImpl::Copy () const
626 { 719 {
627 return new ParetoVariableImpl(*this); 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 {
738 }
642 739
643 //----------------------------------------------------------------------------- 740 // -----------------------------------------------------------------------------
644 //----------------------------------------------------------------------------- 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 v alues,
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 {
806 }
807 WeibullVariableImpl::WeibullVariableImpl (double m, double s)
808 : m_mean (m),
809 m_alpha (s),
810 m_bound (0)
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 }
708 826
709 double WeibullVariableImpl::GetValue() 827 double WeibullVariableImpl::GetValue ()
710 { 828 {
711 if(!m_generator) 829 if (!m_generator)
712 { 830 {
713 m_generator = new RngStream(); 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 {
839 return r;
840 }
841 // otherwise, try again
721 } 842 }
722 } 843 }
723 844
724 RandomVariableBase* WeibullVariableImpl::Copy() const 845 RandomVariableBase* WeibullVariableImpl::Copy () const
725 { 846 {
726 return new WeibullVariableImpl(*this); 847 return new WeibullVariableImpl (*this);
727 } 848 }
728 849
729 WeibullVariable::WeibullVariable() 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 {
865 }
741 866
742 //----------------------------------------------------------------------------- 867 // -----------------------------------------------------------------------------
743 //----------------------------------------------------------------------------- 868 // -----------------------------------------------------------------------------
744 // NormalVariableImpl methods 869 // NormalVariableImpl methods
745 870
746 class NormalVariableImpl : public RandomVariableBase { // Normally Distributed r andom var 871 class NormalVariableImpl : public RandomVariableBase // Normally Distributed r andom 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 me an. 886 * \param b Bound. The NormalVariableImpl is bounded within +-bound of the me an.
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
776 private: 902 private:
777 double m_mean; // Mean value of RV 903 double m_mean; // Mean value of RV
778 double m_variance; // Mean value of RV 904 double m_variance; // Mean value of RV
779 double m_bound; // Bound on value's difference from the mean (absolute val ue) 905 double m_bound; // Bound on value's difference from the mean (absolute val ue)
780 bool m_nextValid; // True if next valid 906 bool m_nextValid; // True if next valid
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),
914 m_variance (1.0),
915 m_bound (INFINITE_VALUE),
916 m_nextValid (false)
917 {
918 }
788 919
789 NormalVariableImpl::NormalVariableImpl(double m, double v, double b) 920 NormalVariableImpl::NormalVariableImpl (double m, double v, double b)
790 : m_mean(m), m_variance(v), m_bound(b), m_nextValid(false) { } 921 : m_mean (m),
922 m_variance (v),
923 m_bound (b),
924 m_nextValid (false)
925 {
926 }
791 927
792 NormalVariableImpl::NormalVariableImpl(const NormalVariableImpl& c) 928 NormalVariableImpl::NormalVariableImpl (const NormalVariableImpl& c)
793 : RandomVariableBase(c), m_mean(c.m_mean), m_variance(c.m_variance), 929 : RandomVariableBase (c),
794 m_bound(c.m_bound), m_nextValid(false) { } 930 m_mean (c.m_mean),
931 m_variance (c.m_variance),
932 m_bound (c.m_bound),
933 m_nextValid (false)
934 {
935 }
795 936
796 double NormalVariableImpl::GetValue() 937 double NormalVariableImpl::GetValue ()
797 { 938 {
798 if(!m_generator) 939 if (!m_generator)
799 { 940 {
800 m_generator = new RngStream(); 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 {
847 return m_mean; 988 return m_mean;
848 } 989 }
849 990
850 double 991 double
851 NormalVariableImpl::GetVariance (void) const 992 NormalVariableImpl::GetVariance (void) const
852 { 993 {
853 return m_variance; 994 return m_variance;
854 } 995 }
855 996
856 double 997 double
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 {
1014 }
871 1015
872 //----------------------------------------------------------------------------- 1016 // -----------------------------------------------------------------------------
873 //----------------------------------------------------------------------------- 1017 // -----------------------------------------------------------------------------
874 class EmpiricalVariableImpl : public RandomVariableBase { 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)
1064 : value (v),
1065 cdf (c)
1066 {
1067 }
1068 EmpiricalVariableImpl::ValueCDF::ValueCDF (const ValueCDF& c)
1069 : value (c.value),
1070 cdf (c.cdf)
1071 {
1072 }
918 1073
919 //----------------------------------------------------------------------------- 1074 // -----------------------------------------------------------------------------
920 //----------------------------------------------------------------------------- 1075 // -----------------------------------------------------------------------------
921 // EmpiricalVariableImpl methods 1076 // EmpiricalVariableImpl methods
922 EmpiricalVariableImpl::EmpiricalVariableImpl() 1077 EmpiricalVariableImpl::EmpiricalVariableImpl ()
923 : validated(false) { } 1078 : validated (false)
1079 {
1080 }
924 1081
925 EmpiricalVariableImpl::EmpiricalVariableImpl(const EmpiricalVariableImpl& c) 1082 EmpiricalVariableImpl::EmpiricalVariableImpl (const EmpiricalVariableImpl& c)
926 : RandomVariableBase(c), validated(c.validated), emp(c.emp) { } 1083 : RandomVariableBase (c),
1084 validated (c.validated),
1085 emp (c.emp)
1086 {
1087 }
927 1088
928 EmpiricalVariableImpl::~EmpiricalVariableImpl() { } 1089 EmpiricalVariableImpl::~EmpiricalVariableImpl ()
1090 {
1091 }
929 1092
930 double EmpiricalVariableImpl::GetValue() 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;
970 } 1133 }
971 else 1134 else
972 { 1135 {
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:
1201 IntEmpiricalVariableImpl ();
1035 1202
1036 IntEmpiricalVariableImpl(); 1203 virtual RandomVariableBase* Copy (void) const;
1037 ··
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
1050 uint32_t IntEmpiricalVariableImpl::GetInteger()
1051 { 1214 {
1052 return (uint32_t)GetValue();
1053 } 1215 }
1054 1216
1055 RandomVariableBase* IntEmpiricalVariableImpl::Copy() const 1217 uint32_t IntEmpiricalVariableImpl::GetInteger ()
1056 { 1218 {
1057 return new IntEmpiricalVariableImpl(*this); 1219 return (uint32_t)GetValue ();
1058 } 1220 }
1059 1221
1060 double IntEmpiricalVariableImpl::Interpolate(double c1, double c2, 1222 RandomVariableBase* IntEmpiricalVariableImpl::Copy () const
1061 double v1, double v2, double r) 1223 {
1062 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2) 1224 return new IntEmpiricalVariableImpl (*this);
1063 return ceil(v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
1064 } 1225 }
1065 1226
1066 IntEmpiricalVariable::IntEmpiricalVariable() 1227 double IntEmpiricalVariableImpl::Interpolate (double c1, double c2,
1228 double v1, double v2, double r)
1229 { // Interpolate random value in range [v1..v2) based on [c1 .. r .. c2)
1230 return ceil (v1 + ((v2 - v1) / (c2 - c1)) * (r - c1));
1231 }
1232
1233 IntEmpiricalVariable::IntEmpiricalVariable ()
1067 : EmpiricalVariable (IntEmpiricalVariableImpl ()) 1234 : EmpiricalVariable (IntEmpiricalVariableImpl ())
1068 {} 1235 {
1236 }
1069 1237
1070 //----------------------------------------------------------------------------- 1238 // -----------------------------------------------------------------------------
1071 //----------------------------------------------------------------------------- 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 ··
1109 double DeterministicVariableImpl::GetValue()
1110 { 1278 {
1111 if (next == count) 1279 }
1280
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 {
1298 }
1126 1299
1127 //----------------------------------------------------------------------------- 1300 // -----------------------------------------------------------------------------
1128 //----------------------------------------------------------------------------- 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 */
1136 LogNormalVariableImpl (double mu, double sigma); 1310 LogNormalVariableImpl (double mu, double sigma);
1137 1311
1138 /** 1312 /**
1139 * \return A random value from this distribution 1313 * \return A random value from this distribution
1140 */ 1314 */
1141 virtual double GetValue (); 1315 virtual double GetValue ();
1142 virtual RandomVariableBase* Copy(void) const; 1316 virtual RandomVariableBase* Copy (void) const;
1143 1317
1144 private: 1318 private:
1145 double m_mu; 1319 double m_mu;
1146 double m_sigma; 1320 double m_sigma;
1147 }; 1321 };
1148 1322
1149 1323
1150 RandomVariableBase* LogNormalVariableImpl::Copy () const 1324 RandomVariableBase* LogNormalVariableImpl::Copy () const
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 {
1197 /* choose x,y in uniform square (-1,-1) to (+1,+1) */ 1372 /* choose x,y in uniform square (-1,-1) to (+1,+1) */
1198 1373
1199 u = -1 + 2 * m_generator->RandU01 (); 1374 u = -1 + 2 * m_generator->RandU01 ();
1200 v = -1 + 2 * m_generator->RandU01 (); 1375 v = -1 + 2 * m_generator->RandU01 ();
1201 1376
1202 /* see if it is in the unit circle */ 1377 /* see if it is in the unit circle */
1203 r2 = u * u + v * v; 1378 r2 = u * u + v * v;
1204 } 1379 }
1205 while (r2 > 1.0 || r2 == 0); 1380 while (r2 > 1.0 || r2 == 0);
1206 1381
1207 normal = u * sqrt (-2.0 * log (r2) / r2); 1382 normal = u * sqrt (-2.0 * log (r2) / r2);
1208 1383
1209 z = exp (m_sigma * normal + m_mu); 1384 z = exp (m_sigma * normal + m_mu);
1210 1385
1211 return z; 1386 return z;
1212 } 1387 }
1213 1388
1214 LogNormalVariable::LogNormalVariable (double mu, double sigma) 1389 LogNormalVariable::LogNormalVariable (double mu, double sigma)
1215 : RandomVariable (LogNormalVariableImpl (mu, sigma)) 1390 : RandomVariable (LogNormalVariableImpl (mu, sigma))
1216 {} 1391 {
1392 }
1217 1393
1218 //----------------------------------------------------------------------------- 1394 // -----------------------------------------------------------------------------
1219 //----------------------------------------------------------------------------- 1395 // -----------------------------------------------------------------------------
1220 // GammaVariableImpl 1396 // GammaVariableImpl
1221 class GammaVariableImpl : public RandomVariableBase 1397 class GammaVariableImpl : public RandomVariableBase
1222 { 1398 {
1223 public: 1399 public:
1224 /** 1400 /**
1225 * \param alpha alpha parameter of the gamma distribution 1401 * \param alpha alpha parameter of the gamma distribution
1226 * \param beta beta parameter of the gamma distribution 1402 * \param beta beta parameter of the gamma distribution
1227 */ 1403 */
1228 GammaVariableImpl (double alpha, double beta); 1404 GammaVariableImpl (double alpha, double beta);
1229 1405
1230 /** 1406 /**
1231 * \return A random value from this distribution 1407 * \return A random value from this distribution
1232 */ 1408 */
1233 virtual double GetValue (); 1409 virtual double GetValue ();
1234 1410
1235 /** 1411 /**
1236 * \return A random value from the gamma distribution with parameters alpha 1412 * \return A random value from the gamma distribution with parameters alpha
1237 * and beta 1413 * and beta
1238 */ 1414 */
1239 double GetValue(double alpha, double beta); 1415 double GetValue (double alpha, double beta);
1240 1416
1241 virtual RandomVariableBase* Copy(void) const; 1417 virtual RandomVariableBase* Copy (void) const;
1242 1418
1243 private: 1419 private:
1244 double m_alpha; 1420 double m_alpha;
1245 double m_beta; 1421 double m_beta;
1246 NormalVariable m_normal; 1422 NormalVariable m_normal;
1247 }; 1423 };
1248 1424
1249 1425
1250 RandomVariableBase* GammaVariableImpl::Copy () const 1426 RandomVariableBase* GammaVariableImpl::Copy () const
1251 { 1427 {
1252 return new GammaVariableImpl (m_alpha, m_beta); 1428 return new GammaVariableImpl (m_alpha, m_beta);
1253 } 1429 }
1254 1430
1255 GammaVariableImpl::GammaVariableImpl (double alpha, double beta) 1431 GammaVariableImpl::GammaVariableImpl (double alpha, double beta)
1256 : m_alpha(alpha), m_beta(beta) 1432 : m_alpha (alpha),
1433 m_beta (beta)
1257 { 1434 {
1258 } 1435 }
1259 1436
1260 double 1437 double
1261 GammaVariableImpl::GetValue () 1438 GammaVariableImpl::GetValue ()
1262 { 1439 {
1263 return GetValue(m_alpha, m_beta); 1440 return GetValue (m_alpha, m_beta);
1264 } 1441 }
1265 1442
1266 /* 1443 /*
1267 The code for the following generator functions was adapted from ns-2 1444 The code for the following generator functions was adapted from ns-2
1268 tools/ranvar.cc 1445 tools/ranvar.cc
1269 1446
1270 Originally the algorithm was devised by Marsaglia in 2000: 1447 Originally the algorithm was devised by Marsaglia in 2000:
1271 G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables 1448 G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
1272 ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000 1449 ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
1273 1450
1274 The Gamma distribution density function has the form 1451 The Gamma distribution density function has the form
1275 1452
1276 » x^(alpha-1) * exp(-x/beta) 1453 x^(alpha-1) * exp(-x/beta)
1277 » p(x; alpha, beta) = ---------------------------- 1454 p(x; alpha, beta) = ----------------------------
1278 beta^alpha * Gamma(alpha) 1455 beta^alpha * Gamma(alpha)
1279 1456
1280 for x > 0. 1457 for x > 0.
1281 */ 1458 */
1282 double 1459 double
1283 GammaVariableImpl::GetValue (double alpha, double beta) 1460 GammaVariableImpl::GetValue (double alpha, double beta)
1284 { 1461 {
1285 if(!m_generator) 1462 if (!m_generator)
1286 { 1463 {
1287 m_generator = new RngStream(); 1464 m_generator = new RngStream ();
1288 } 1465 }
1289 1466
1290 if (alpha < 1) 1467 if (alpha < 1)
1291 { 1468 {
1292 double u = m_generator->RandU01 (); 1469 double u = m_generator->RandU01 ();
1293 return GetValue(1.0 + alpha, beta) * pow (u, 1.0 / alpha); 1470 return GetValue (1.0 + alpha, beta) * pow (u, 1.0 / alpha);
1294 } 1471 }
1295 » 1472
1296 double x, v, u; 1473 double x, v, u;
1297 double d = alpha - 1.0 / 3.0; 1474 double d = alpha - 1.0 / 3.0;
1298 double c = (1.0 / 3.0) / sqrt (d); 1475 double c = (1.0 / 3.0) / sqrt (d);
1299 1476
1300 while (1) 1477 while (1)
1301 { 1478 {
1302 do 1479 do
1303 { 1480 {
1304 x = m_normal.GetValue (); 1481 x = m_normal.GetValue ();
1305 v = 1.0 + c * x; 1482 v = 1.0 + c * x;
1306 } while (v <= 0); 1483 }
1484 while (v <= 0);
1307 1485
1308 v = v * v * v; 1486 v = v * v * v;
1309 u = m_generator->RandU01 (); 1487 u = m_generator->RandU01 ();
1310 if (u < 1 - 0.0331 * x * x * x * x) 1488 if (u < 1 - 0.0331 * x * x * x * x)
1311 { 1489 {
1312 break; 1490 break;
1313 } 1491 }
1314 if (log (u) < 0.5 * x * x + d * (1 - v + log (v))) 1492 if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
1315 { 1493 {
1316 break; 1494 break;
1317 } 1495 }
1318 } 1496 }
1319 1497
1320 return beta * d * v; 1498 return beta * d * v;
1321 } 1499 }
1322 1500
1323 GammaVariable::GammaVariable () 1501 GammaVariable::GammaVariable ()
1324 : RandomVariable (GammaVariableImpl (1.0, 1.0)) 1502 : RandomVariable (GammaVariableImpl (1.0, 1.0))
1325 { 1503 {
1326 } 1504 }
1327 1505
1328 GammaVariable::GammaVariable (double alpha, double beta) 1506 GammaVariable::GammaVariable (double alpha, double beta)
1329 : RandomVariable (GammaVariableImpl (alpha, beta)) 1507 : RandomVariable (GammaVariableImpl (alpha, beta))
1330 { 1508 {
1331 } 1509 }
1332 1510
1333 double GammaVariable::GetValue(void) const 1511 double GammaVariable::GetValue (void) const
1334 { 1512 {
1335 return this->RandomVariable::GetValue (); 1513 return this->RandomVariable::GetValue ();
1336 } 1514 }
1337 1515
1338 double GammaVariable::GetValue(double alpha, double beta) const 1516 double GammaVariable::GetValue (double alpha, double beta) const
1339 { 1517 {
1340 return ((GammaVariableImpl*)Peek())->GetValue(alpha, beta); 1518 return ((GammaVariableImpl*)Peek ())->GetValue (alpha, beta);
1341 } 1519 }
1342 1520
1343 //----------------------------------------------------------------------------- 1521 // -----------------------------------------------------------------------------
1344 //----------------------------------------------------------------------------- 1522 // -----------------------------------------------------------------------------
1345 // ErlangVariableImpl 1523 // ErlangVariableImpl
1346 1524
1347 class ErlangVariableImpl : public RandomVariableBase 1525 class ErlangVariableImpl : public RandomVariableBase
1348 { 1526 {
1349 public: 1527 public:
1350 /** 1528 /**
1351 * \param k k parameter of the Erlang distribution 1529 * \param k k parameter of the Erlang distribution
1352 * \param lambda lambda parameter of the Erlang distribution 1530 * \param lambda lambda parameter of the Erlang distribution
1353 */ 1531 */
1354 ErlangVariableImpl (unsigned int k, double lambda); 1532 ErlangVariableImpl (unsigned int k, double lambda);
1355 1533
1356 /** 1534 /**
1357 * \return A random value from this distribution 1535 * \return A random value from this distribution
1358 */ 1536 */
1359 virtual double GetValue (); 1537 virtual double GetValue ();
1360 1538
1361 /** 1539 /**
1362 * \return A random value from the Erlang distribution with parameters k and 1540 * \return A random value from the Erlang distribution with parameters k and
1363 * lambda. 1541 * lambda.
1364 */ 1542 */
1365 double GetValue(unsigned int k, double lambda); 1543 double GetValue (unsigned int k, double lambda);
1366 1544
1367 virtual RandomVariableBase* Copy(void) const; 1545 virtual RandomVariableBase* Copy (void) const;
1368 1546
1369 private: 1547 private:
1370 unsigned int m_k; 1548 unsigned int m_k;
1371 double m_lambda; 1549 double m_lambda;
1372 }; 1550 };
1373 1551
1374 1552
1375 RandomVariableBase* ErlangVariableImpl::Copy () const 1553 RandomVariableBase* ErlangVariableImpl::Copy () const
1376 { 1554 {
1377 return new ErlangVariableImpl (m_k, m_lambda); 1555 return new ErlangVariableImpl (m_k, m_lambda);
1378 } 1556 }
1379 1557
1380 ErlangVariableImpl::ErlangVariableImpl (unsigned int k, double lambda) 1558 ErlangVariableImpl::ErlangVariableImpl (unsigned int k, double lambda)
1381 : m_k(k), m_lambda(lambda) 1559 : m_k (k),
1560 m_lambda (lambda)
1382 { 1561 {
1383 } 1562 }
1384 1563
1385 double 1564 double
1386 ErlangVariableImpl::GetValue () 1565 ErlangVariableImpl::GetValue ()
1387 { 1566 {
1388 return GetValue(m_k, m_lambda); 1567 return GetValue (m_k, m_lambda);
1389 } 1568 }
1390 1569
1391 /* 1570 /*
1392 The code for the following generator functions was adapted from ns-2 1571 The code for the following generator functions was adapted from ns-2
1393 tools/ranvar.cc 1572 tools/ranvar.cc
1394 1573
1395 The Erlang distribution density function has the form 1574 The Erlang distribution density function has the form
1396 1575
1397 » x^(k-1) * exp(-x/lambda) 1576 x^(k-1) * exp(-x/lambda)
1398 » p(x; k, lambda) = --------------------------- 1577 p(x; k, lambda) = ---------------------------
1399 lambda^k * (k-1)! 1578 lambda^k * (k-1)!
1400 1579
1401 for x > 0. 1580 for x > 0.
1402 */ 1581 */
1403 double 1582 double
1404 ErlangVariableImpl::GetValue (unsigned int k, double lambda) 1583 ErlangVariableImpl::GetValue (unsigned int k, double lambda)
1405 { 1584 {
1406 if(!m_generator) 1585 if (!m_generator)
1407 { 1586 {
1408 m_generator = new RngStream(); 1587 m_generator = new RngStream ();
1409 } 1588 }
1410 1589
1411 ExponentialVariable exponential(lambda); 1590 ExponentialVariable exponential (lambda);
1412 1591
1413 double result = 0; 1592 double result = 0;
1414 for (unsigned int i = 0; i < k; ++i) 1593 for (unsigned int i = 0; i < k; ++i)
1415 { 1594 {
1416 result += exponential.GetValue(); 1595 result += exponential.GetValue ();
1417 } 1596 }
1418 1597
1419 return result; 1598 return result;
1420 } 1599 }
1421 1600
1422 ErlangVariable::ErlangVariable () 1601 ErlangVariable::ErlangVariable ()
1423 : RandomVariable (ErlangVariableImpl (1, 1.0)) 1602 : RandomVariable (ErlangVariableImpl (1, 1.0))
1424 { 1603 {
1425 } 1604 }
1426 1605
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)·
1484 : RandomVariableBase(c), m_min(c.m_min), m_max(c.m_max), m_mode(c.m_mode) { }
1485 1663
1486 double TriangularVariableImpl::GetValue() 1664 TriangularVariableImpl::TriangularVariableImpl (double s, double l, double mean)
1665 : m_min (s),
1666 m_max (l),
1667 m_mode (3.0 * mean - s - l)
1487 { 1668 {
1488 if(!m_generator) 1669 }
1670
1671 TriangularVariableImpl::TriangularVariableImpl (const TriangularVariableImpl& c)
1672 : RandomVariableBase (c),
1673 m_min (c.m_min),
1674 m_max (c.m_max),
1675 m_mode (c.m_mode)
1676 {
1677 }
1678
1679 double TriangularVariableImpl::GetValue ()
1680 {
1681 if (!m_generator)
1489 { 1682 {
1490 m_generator = new RngStream(); 1683 m_generator = new RngStream ();
1491 } 1684 }
1492 double u = m_generator->RandU01(); 1685 double u = m_generator->RandU01 ();
1493 if(u <= (m_mode - m_min) / (m_max - m_min) ) 1686 if (u <= (m_mode - m_min) / (m_max - m_min) )
1494 { 1687 {
1495 return m_min + sqrt(u * (m_max - m_min) * (m_mode - m_min) ); 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 {
1708 }
1514 1709
1515 //----------------------------------------------------------------------------- 1710 // -----------------------------------------------------------------------------
1516 //----------------------------------------------------------------------------- 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 */
1524 ZipfVariableImpl (long n, double alpha); 1720 ZipfVariableImpl (long n, double alpha);
1525 1721
1526 /** 1722 /**
1527 * \A zipf variable with N=1 and alpha=0 1723 * \A zipf variable with N=1 and alpha=0
1528 */ 1724 */
1529 ZipfVariableImpl (); 1725 ZipfVariableImpl ();
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),
1755 m_alpha (alpha),
1756 m_c (0)
1557 { 1757 {
1558 //calculate the normalization constant c 1758 // calculate the normalization constant c
1559 for(int i=1;i<=n;i++) 1759 for (int i = 1; i <= n; i++)
1560 { 1760 {
1561 m_c+=(1.0/pow((double)i,alpha)); 1761 m_c += (1.0 / pow ((double)i,alpha));
1562 } 1762 }
1563 m_c=1.0/m_c; 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 {
1796 }
1595 1797
1596 1798
1597 std::ostream &operator << (std::ostream &os, const RandomVariable &var) 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, 1 986.
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 {
1603 os << "Constant:" << constant->GetValue (); 1892 os << "Constant:" << constant->GetValue ();
1604 return os; 1893 return os;
1605 } 1894 }
1606 UniformVariableImpl *uniform = dynamic_cast<UniformVariableImpl *> (base); 1895 UniformVariableImpl *uniform = dynamic_cast<UniformVariableImpl *> (base);
1607 if (uniform != 0) 1896 if (uniform != 0)
1608 { 1897 {
1609 os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax (); 1898 os << "Uniform:" << uniform->GetMin () << ":" << uniform->GetMax ();
1610 return os; 1899 return os;
1611 } 1900 }
1612 NormalVariableImpl *normal = dynamic_cast<NormalVariableImpl *> (base); 1901 NormalVariableImpl *normal = dynamic_cast<NormalVariableImpl *> (base);
1613 if (normal != 0) 1902 if (normal != 0)
1614 { 1903 {
1615 os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance (); 1904 os << "Normal:" << normal->GetMean () << ":" << normal->GetVariance ();
1616 double bound = normal->GetBound (); 1905 double bound = normal->GetBound ();
1617 if (bound != NormalVariableImpl::INFINITE_VALUE) 1906 if (bound != NormalVariableImpl::INFINITE_VALUE)
1618 { 1907 {
1619 os << ":" << bound; 1908 os << ":" << bound;
1620 } 1909 }
1621 return os; 1910 return os;
1622 } 1911 }
1623 // XXX: support other distributions 1912 // XXX: support other distributions
1624 os.setstate (std::ios_base::badbit); 1913 os.setstate (std::ios_base::badbit);
1625 return os; 1914 return os;
1626 } 1915 }
1627 std::istream &operator >> (std::istream &is, RandomVariable &var) 1916 std::istream & operator >> (std::istream &is, RandomVariable &var)
1628 { 1917 {
1629 std::string value; 1918 std::string value;
1630 is >> value; 1919 is >> value;
1631 std::string::size_type tmp; 1920 std::string::size_type tmp;
1632 tmp = value.find (":"); 1921 tmp = value.find (":");
1633 if (tmp == std::string::npos) 1922 if (tmp == std::string::npos)
1634 { 1923 {
1635 is.setstate (std::ios_base::badbit); 1924 is.setstate (std::ios_base::badbit);
1636 return is; 1925 return is;
1637 } 1926 }
(...skipping 70 matching lines...) Expand 10 before | Expand all | Expand 10 after
1708 NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << t ype); 1997 NS_FATAL_ERROR ("RandomVariable deserialization not implemented for " << t ype);
1709 // XXX: support other distributions. 1998 // XXX: support other distributions.
1710 } 1999 }
1711 return is; 2000 return is;
1712 } 2001 }
1713 2002
1714 class BasicRandomNumberTestCase : public TestCase 2003 class BasicRandomNumberTestCase : public TestCase
1715 { 2004 {
1716 public: 2005 public:
1717 BasicRandomNumberTestCase (); 2006 BasicRandomNumberTestCase ();
1718 virtual ~BasicRandomNumberTestCase () {} 2007 virtual ~BasicRandomNumberTestCase ()
2008 {
2009 }
1719 2010
1720 private: 2011 private:
1721 virtual bool DoRun (void); 2012 virtual bool DoRun (void);
1722 }; 2013 };
1723 2014
1724 BasicRandomNumberTestCase::BasicRandomNumberTestCase () 2015 BasicRandomNumberTestCase::BasicRandomNumberTestCase ()
1725 : TestCase ("Check basic random number operation") 2016 : TestCase ("Check basic random number operation")
1726 { 2017 {
1727 } 2018 }
1728 2019
1729 bool 2020 bool
1730 BasicRandomNumberTestCase::DoRun (void) 2021 BasicRandomNumberTestCase::DoRun (void)
1731 { 2022 {
1732 const double desiredMean = 1.0; 2023 const double desiredMean = 1.0;
1733 const double desiredStdDev = 1.0; 2024 const double desiredStdDev = 1.0;
1734 2025
1735 double tmp = log (1 + (desiredStdDev / desiredMean) * (desiredStdDev / desired Mean)); 2026 double tmp = log (1 + (desiredStdDev / desiredMean) * (desiredStdDev / desired Mean));
1736 double sigma = sqrt (tmp); 2027 double sigma = sqrt (tmp);
1737 double mu = log (desiredMean) - 0.5 * tmp; 2028 double mu = log (desiredMean) - 0.5 * tmp;
1738 2029
1739 // 2030 //
1740 // Test a custom lognormal instance to see if its moments have any relation 2031 // Test a custom lognormal instance to see if its moments have any relation
1741 // expected reality. 2032 // expected reality.
1742 // 2033 //
1743 LogNormalVariable lognormal (mu, sigma); 2034 LogNormalVariable lognormal (mu, sigma);
1744 vector<double> samples; 2035 vector<double> samples;
1745 const int NSAMPLES = 10000; 2036 const int NSAMPLES = 10000;
1746 double sum = 0; 2037 double sum = 0;
1747 2038
1748 // 2039 //
1749 // Get and store a bunch of samples. As we go along sum them and then find 2040 // Get and store a bunch of samples. As we go along sum them and then find
1750 // the mean value of the samples. 2041 // the mean value of the samples.
1751 // 2042 //
1752 for (int n = NSAMPLES; n; --n) 2043 for (int n = NSAMPLES; n; --n)
1753 { 2044 {
1754 double value = lognormal.GetValue (); 2045 double value = lognormal.GetValue ();
1755 sum += value; 2046 sum += value;
1756 samples.push_back (value); 2047 samples.push_back (value);
1757 } 2048 }
1758 double obtainedMean = sum / NSAMPLES; 2049 double obtainedMean = sum / NSAMPLES;
1759 NS_TEST_EXPECT_MSG_EQ_TOL (obtainedMean, desiredMean, 0.1, "Got unexpected mea n value from LogNormalVariable"); 2050 NS_TEST_EXPECT_MSG_EQ_TOL (obtainedMean, desiredMean, 0.1, "Got unexpected mea n value from LogNormalVariable");
1760 2051
1761 // 2052 //
1762 // Wander back through the saved stamples and find their standard deviation 2053 // Wander back through the saved stamples and find their standard deviation
1763 // 2054 //
1764 sum = 0; 2055 sum = 0;
1765 for (vector<double>::iterator iter = samples.begin (); iter != samples.end (); iter++) 2056 for (vector<double>::iterator iter = samples.begin (); iter != samples.end (); iter++)
1766 { 2057 {
1767 double tmp = (*iter - obtainedMean); 2058 double tmp = (*iter - obtainedMean);
1768 sum += tmp * tmp; 2059 sum += tmp * tmp;
1769 } 2060 }
1770 double obtainedStdDev = sqrt (sum / (NSAMPLES - 1)); 2061 double obtainedStdDev = sqrt (sum / (NSAMPLES - 1));
1771 NS_TEST_EXPECT_MSG_EQ_TOL (obtainedStdDev, desiredStdDev, 0.1, "Got unexpected standard deviation from LogNormalVariable"); 2062 NS_TEST_EXPECT_MSG_EQ_TOL (obtainedStdDev, desiredStdDev, 0.1, "Got unexpected standard deviation from LogNormalVariable");
1772 2063
1773 return GetErrorStatus (); 2064 return GetErrorStatus ();
1774 } 2065 }
1775 2066
1776 class RandomNumberSerializationTestCase : public TestCase 2067 class RandomNumberSerializationTestCase : public TestCase
1777 { 2068 {
1778 public: 2069 public:
1779 RandomNumberSerializationTestCase (); 2070 RandomNumberSerializationTestCase ();
1780 virtual ~RandomNumberSerializationTestCase () {} 2071 virtual ~RandomNumberSerializationTestCase ()
2072 {
2073 }
1781 2074
1782 private: 2075 private:
1783 virtual bool DoRun (void); 2076 virtual bool DoRun (void);
1784 }; 2077 };
1785 2078
1786 RandomNumberSerializationTestCase::RandomNumberSerializationTestCase () 2079 RandomNumberSerializationTestCase::RandomNumberSerializationTestCase ()
1787 : TestCase ("Check basic random number operation") 2080 : TestCase ("Check basic random number operation")
1788 { 2081 {
1789 } 2082 }
1790 2083
1791 bool 2084 bool
1792 RandomNumberSerializationTestCase::DoRun (void) 2085 RandomNumberSerializationTestCase::DoRun (void)
1793 { 2086 {
1794 RandomVariableValue val; 2087 RandomVariableValue val;
1795 val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ()); 2088 val.DeserializeFromString ("Uniform:0.1:0.2", MakeRandomVariableChecker ());
1796 RandomVariable rng = val.Get (); 2089 RandomVariable rng = val.Get ();
1797 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Uniform:0.1:0.2", 2090 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Uniform:0.1:0.2",
1798 "Deserialize and Serialize \"Uniform:0.1:0.2\" mismatch "); 2091 "Deserialize and Serialize \"Uniform:0.1:0.2\" mismatch ");
1799 2092
1800 val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ()); 2093 val.DeserializeFromString ("Normal:0.1:0.2", MakeRandomVariableChecker ());
1801 rng = val.Get (); 2094 rng = val.Get ();
1802 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Normal:0.1:0.2", 2095 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Normal:0.1:0.2",
1803 "Deserialize and Serialize \"Normal:0.1:0.2\" mismatch" ); 2096 "Deserialize and Serialize \"Normal:0.1:0.2\" mismatch" );
1804 2097
1805 val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker () ); 2098 val.DeserializeFromString ("Normal:0.1:0.2:0.15", MakeRandomVariableChecker () );
1806 rng = val.Get (); 2099 rng = val.Get ();
1807 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Normal:0.1:0.2:0.15", 2100 NS_TEST_ASSERT_MSG_EQ (val.SerializeToString (MakeRandomVariableChecker ()), " Normal:0.1:0.2:0.15",
(...skipping 10 matching lines...) Expand all
1818 2111
1819 BasicRandomNumberTestSuite::BasicRandomNumberTestSuite () 2112 BasicRandomNumberTestSuite::BasicRandomNumberTestSuite ()
1820 : TestSuite ("basic-random-number", BVT) 2113 : TestSuite ("basic-random-number", BVT)
1821 { 2114 {
1822 AddTestCase (new BasicRandomNumberTestCase); 2115 AddTestCase (new BasicRandomNumberTestCase);
1823 AddTestCase (new RandomNumberSerializationTestCase); 2116 AddTestCase (new RandomNumberSerializationTestCase);
1824 } 2117 }
1825 2118
1826 BasicRandomNumberTestSuite BasicRandomNumberTestSuite; 2119 BasicRandomNumberTestSuite BasicRandomNumberTestSuite;
1827 2120
1828 }//namespace ns3 2121 } // namespace ns3
OLDNEW
« no previous file with comments | « src/core/random-variable.h ('k') | no next file » | no next file with comments »

Powered by Google App Engine
RSS Feeds Recent Issues | This issue
This is Rietveld f62528b