[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

random.hxx
1/************************************************************************/
2/* */
3/* Copyright 2008 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_RANDOM_HXX
38#define VIGRA_RANDOM_HXX
39
40#include "mathutil.hxx"
41#include "functortraits.hxx"
42#include "array_vector.hxx"
43
44#include <ctime>
45#include <cstdint>
46
47 // includes to get the current process and thread IDs
48 // to be used for automated seeding
49#ifdef _MSC_VER
50 #include <vigra/windows.h> // for GetCurrentProcessId() and GetCurrentThreadId()
51#endif
52
53#ifdef __linux__
54 #include <unistd.h> // for getpid()
55 #include <sys/syscall.h> // for SYS_gettid
56#endif
57
58#ifdef __APPLE__
59 #include <pthread.h>
60 #include <unistd.h> // for getpid()
61 #include <sys/syscall.h> // SYS_thread_selfid
62 #include <AvailabilityMacros.h> // to check if we are on MacOS 10.6 or later
63#endif
64
65namespace vigra {
66
67enum RandomSeedTag { RandomSeed };
68
69namespace detail {
70
71enum RandomEngineTag { TT800, MT19937 };
72
73
74template<RandomEngineTag EngineTag>
75struct RandomState;
76
77template <RandomEngineTag EngineTag>
78void seed(UInt32 theSeed, RandomState<EngineTag> & engine)
79{
80 engine.state_[0] = theSeed;
81 for(UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
82 {
83 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
84 }
85}
86
87template <class Iterator, RandomEngineTag EngineTag>
88void seed(Iterator init, UInt32 key_length, RandomState<EngineTag> & engine)
89{
90 const UInt32 N = RandomState<EngineTag>::N;
91 int k = static_cast<int>(std::max(N, key_length));
92 UInt32 i = 1, j = 0;
93 Iterator data = init;
94 for (; k; --k)
95 {
96 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
97 + *data + j; /* non linear */
98 ++i; ++j; ++data;
99
100 if (i >= N)
101 {
102 engine.state_[0] = engine.state_[N-1];
103 i=1;
104 }
105 if (j>=key_length)
106 {
107 j=0;
108 data = init;
109 }
110 }
111
112 for (k=N-1; k; --k)
113 {
114 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
115 - i; /* non linear */
116 ++i;
117 if (i>=N)
118 {
119 engine.state_[0] = engine.state_[N-1];
120 i=1;
121 }
122 }
123
124 engine.state_[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
125}
126
127template <RandomEngineTag EngineTag>
128void seed(RandomSeedTag, RandomState<EngineTag> & engine)
129{
130 static UInt32 globalCount = 0;
131 ArrayVector<UInt32> seedData;
132
133 seedData.push_back(static_cast<UInt32>(time(0)));
134 seedData.push_back(static_cast<UInt32>(clock()));
135 seedData.push_back(++globalCount);
136
137 const uintptr_t ptr = reinterpret_cast<uintptr_t>(&engine);
138 seedData.push_back(static_cast<UInt32>((ptr & 0xffffffff)));
139 static const UInt32 shift = sizeof(ptr) > 4 ? 32 : 16;
140 seedData.push_back(static_cast<UInt32>((ptr >> shift)));
141
142#ifdef _MSC_VER
143 seedData.push_back(static_cast<UInt32>(GetCurrentProcessId()));
144 seedData.push_back(static_cast<UInt32>(GetCurrentThreadId()));
145#endif
146
147#ifdef __linux__
148 seedData.push_back(static_cast<UInt32>(getpid()));
149# ifdef SYS_gettid
150 seedData.push_back(static_cast<UInt32>(syscall(SYS_gettid)));
151# endif
152#endif
153
154#ifdef __APPLE__
155 seedData.push_back(static_cast<UInt32>(getpid()));
156 #if __MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_12
157 uint64_t tid64;
158 pthread_threadid_np(NULL, &tid64);
159 seedData.push_back(static_cast<UInt32>(tid64));
160 #elif defined(SYS_thread_selfid) && (__MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_6)
161 // SYS_thread_selfid was introduced in MacOS 10.6
162 seedData.push_back(static_cast<UInt32>(syscall(SYS_thread_selfid)));
163 #elif defined(SYS_gettid)
164 seedData.push_back(static_cast<UInt32>(syscall(SYS_gettid)));
165 #endif
166#endif
167
168 seed(seedData.begin(), seedData.size(), engine);
169}
170
171 /* Tempered twister TT800 by M. Matsumoto */
172template<>
173struct RandomState<TT800>
174{
175 static const UInt32 N = 25, M = 7;
176
177 mutable UInt32 state_[N];
178 mutable UInt32 current_;
179
180 RandomState()
181 : current_(0)
182 {
183 UInt32 seeds[N] = {
184 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
185 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
186 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
187 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
188 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
189 };
190
191 for(UInt32 i=0; i<N; ++i)
192 state_[i] = seeds[i];
193 }
194
195 protected:
196
197 UInt32 get() const
198 {
199 if(current_ == N)
200 generateNumbers<void>();
201
202 UInt32 y = state_[current_++];
203 y ^= (y << 7) & 0x2b5b2500;
204 y ^= (y << 15) & 0xdb8b0000;
205 return y ^ (y >> 16);
206 }
207
208 template <class DUMMY>
209 void generateNumbers() const;
210
211 void seedImpl(RandomSeedTag)
212 {
213 seed(RandomSeed, *this);
214 }
215
216 void seedImpl(UInt32 theSeed)
217 {
218 seed(theSeed, *this);
219 }
220
221 template<class Iterator>
222 void seedImpl(Iterator init, UInt32 length)
223 {
224 seed(init, length, *this);
225 }
226};
227
228template <class DUMMY>
229void RandomState<TT800>::generateNumbers() const
230{
231 UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
232
233 for(UInt32 i=0; i<N-M; ++i)
234 {
235 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
236 }
237 for (UInt32 i=N-M; i<N; ++i)
238 {
239 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
240 }
241 current_ = 0;
242}
243
244 /* Mersenne twister MT19937 by M. Matsumoto */
245template<>
246struct RandomState<MT19937>
247{
248 static const UInt32 N = 624, M = 397;
249
250 mutable UInt32 state_[N];
251 mutable UInt32 current_;
252
253 RandomState()
254 : current_(0)
255 {
256 seed(19650218U, *this);
257 }
258
259 protected:
260
261 UInt32 get() const
262 {
263 if(current_ == N)
264 generateNumbers<void>();
265
266 UInt32 x = state_[current_++];
267 x ^= (x >> 11);
268 x ^= (x << 7) & 0x9D2C5680U;
269 x ^= (x << 15) & 0xEFC60000U;
270 return x ^ (x >> 18);
271 }
272
273 template <class DUMMY>
274 void generateNumbers() const;
275
276 static UInt32 twiddle(UInt32 u, UInt32 v)
277 {
278 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
279 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
280 }
281
282 void seedImpl(RandomSeedTag)
283 {
284 seed(RandomSeed, *this);
285 generateNumbers<void>();
286 }
287
288 void seedImpl(UInt32 theSeed)
289 {
290 seed(theSeed, *this);
291 generateNumbers<void>();
292 }
293
294 template<class Iterator>
295 void seedImpl(Iterator init, UInt32 length)
296 {
297 seed(19650218U, *this);
298 seed(init, length, *this);
299 generateNumbers<void>();
300 }
301};
302
303template <class DUMMY>
304void RandomState<MT19937>::generateNumbers() const
305{
306 for (unsigned int i = 0; i < (N - M); ++i)
307 {
308 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
309 }
310 for (unsigned int i = N - M; i < (N - 1); ++i)
311 {
312 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
313 }
314 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
315 current_ = 0;
316}
317
318} // namespace detail
319
320
321/** \addtogroup RandomNumberGeneration Random Number Generation
322
323 High-quality random number generators and related functors.
324*/
325//@{
326
327/** Generic random number generator.
328
329 The actual generator is passed in the template argument <tt>Engine</tt>. Two generators
330 are currently available:
331 <ul>
332 <li> <tt>RandomMT19937</tt>: The state-of-the-art <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">Mersenne Twister</a> with a state length of 2<sup>19937</sup> and very high statistical quality.
333 <li> <tt>RandomTT800</tt>: (default) The Tempered Twister, a simpler predecessor of the Mersenne Twister with period length 2<sup>800</sup>.
334 </ul>
335
336 Both generators have been designed by <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html">Makoto Matsumoto</a>.
337
338 <b>Traits defined:</b>
339
340 \verbatim FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer \endverbatim
341 is true (<tt>VigraTrueType</tt>).
342*/
343template <class Engine = detail::RandomState<detail::MT19937> >
345: public Engine
346{
347 mutable double normalCached_;
348 mutable bool normalCachedValid_;
349
350 public:
351
352 /** Create a new random generator object with standard seed.
353
354 Due to standard seeding, the random numbers generated will always be the same.
355 This is useful for debugging.
356 */
358 : normalCached_(0.0),
359 normalCachedValid_(false)
360 {}
361
362 /** Create a new random generator object with a random seed.
363
364 The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
365 values.
366
367 <b>Usage:</b>
368 \code
369 RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed);
370 \endcode
371 */
373 : normalCached_(0.0),
374 normalCachedValid_(false)
375 {
376 this->seedImpl(RandomSeed);
377 }
378
379 /** Create a new random generator object from the given seed.
380
381 The same seed will always produce identical random sequences.
382 If \a ignoreSeed is <tt>true</tt>, the given seed is ignored,
383 and the generator is seeded randomly (as if it was constructed
384 with <tt>RandomNumberGenerator<>(RandomSeed)</tt>). This allows
385 you to switch between random and deterministic seeding at
386 run-time.
387 */
388 RandomNumberGenerator(UInt32 theSeed, bool ignoreSeed=false)
389 : normalCached_(0.0),
390 normalCachedValid_(false)
391 {
392 if(ignoreSeed)
393 this->seedImpl(RandomSeed);
394 else
395 this->seedImpl(theSeed);
396 }
397
398 /** Create a new random generator object from the given seed sequence.
399
400 Longer seed sequences lead to better initialization in the sense that the generator's
401 state space is covered much better than is possible with 32-bit seeds alone.
402 */
403 template<class Iterator>
404 RandomNumberGenerator(Iterator init, UInt32 length)
405 : normalCached_(0.0),
406 normalCachedValid_(false)
407 {
408 this->seedImpl(init, length);
409 }
410
411
412 /** Re-initialize the random generator object with a random seed.
413
414 The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
415 values.
416
417 <b>Usage:</b>
418 \code
419 RandomNumberGenerator<> rnd = ...;
420 ...
421 rnd.seed(RandomSeed);
422 \endcode
423 */
424 void seed(RandomSeedTag)
425 {
426 this->seedImpl(RandomSeed);
427 normalCachedValid_ = false;
428 }
429
430 /** Re-initialize the random generator object from the given seed.
431
432 The same seed will always produce identical random sequences.
433 If \a ignoreSeed is <tt>true</tt>, the given seed is ignored,
434 and the generator is seeded randomly (as if <tt>seed(RandomSeed)</tt>
435 was called). This allows you to switch between random and deterministic
436 seeding at run-time.
437 */
438 void seed(UInt32 theSeed, bool ignoreSeed=false)
439 {
440 if(ignoreSeed)
441 this->seedImpl(RandomSeed);
442 else
443 this->seedImpl(theSeed);
444 normalCachedValid_ = false;
445 }
446
447 /** Re-initialize the random generator object from the given seed sequence.
448
449 Longer seed sequences lead to better initialization in the sense that the generator's
450 state space is covered much better than is possible with 32-bit seeds alone.
451 */
452 template<class Iterator>
453 void seed(Iterator init, UInt32 length)
454 {
455 this->seedImpl(init, length);
456 normalCachedValid_ = false;
457 }
458
459 /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
460
461 That is, 0 &lt;= i &lt; 2<sup>32</sup>.
462 */
464 {
465 return this->get();
466 }
467
468 /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
469
470 That is, 0 &lt;= i &lt; 2<sup>32</sup>.
471 */
473 {
474 return this->get();
475 }
476
477
478#if 0 // difficult implementation necessary if low bits are not sufficiently random
479 // in [0,beyond)
480 UInt32 uniformInt(UInt32 beyond) const
481 {
482 if(beyond < 2)
483 return 0;
484
485 UInt32 factor = factorForUniformInt(beyond);
486 UInt32 res = this->get() / factor;
487
488 // Use rejection method to avoid quantization bias.
489 // On average, we will need two raw random numbers to generate one.
490 while(res >= beyond)
491 res = this->get() / factor;
492 return res;
493 }
494#endif /* #if 0 */
495
496 /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>).
497
498 That is, 0 &lt;= i &lt; <tt>beyond</tt>.
499 */
501 {
502 // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is
503 // the case for TT800 and MT19937
504 if(beyond < 2)
505 return 0;
506
507 UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
508 UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
509 UInt32 res = this->get();
510
511 // Use rejection method to avoid quantization bias.
512 // We will need two raw random numbers in amortized worst case.
513 while(res > lastSafeValue)
514 res = this->get();
515 return res % beyond;
516 }
517
518 /** Return a uniformly distributed double-precision random number in [0.0, 1.0).
519
520 That is, 0.0 &lt;= i &lt; 1.0. All 53-bit bits of the mantissa are random (two 32-bit integers are used to
521 create this number).
522 */
523 double uniform53() const
524 {
525 // make full use of the entire 53-bit mantissa of a double, by Isaku Wada
526 return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0);
527 }
528
529 /** Return a uniformly distributed double-precision random number in [0.0, 1.0].
530
531 That is, 0.0 &lt;= i &lt;= 1.0. This number is computed by <tt>uniformInt()</tt> / (2<sup>32</sup> - 1),
532 so it has effectively only 32 random bits.
533 */
534 double uniform() const
535 {
536 return static_cast<double>(this->get()) / 4294967295.0;
537 }
538
539 /** Return a uniformly distributed double-precision random number in [lower, upper].
540
541 That is, <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>. This number is computed
542 from <tt>uniform()</tt>, so it has effectively only 32 random bits.
543 */
544 double uniform(double lower, double upper) const
545 {
546 vigra_precondition(lower < upper,
547 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
548 return uniform() * (upper-lower) + lower;
549 }
550
551 /** Return a standard normal variate (Gaussian) random number.
552
553 Mean is zero, standard deviation is 1.0. It uses the polar form of the
554 Box-Muller transform.
555 */
556 double normal() const;
557
558 /** Return a normal variate (Gaussian) random number with the given mean and standard deviation.
559
560 It uses the polar form of the Box-Muller transform.
561 */
562 double normal(double mean, double stddev) const
563 {
564 vigra_precondition(stddev > 0.0,
565 "RandomNumberGenerator::normal(): standard deviation must be positive.");
566 return normal()*stddev + mean;
567 }
568
569 /** Access the global (program-wide) instance of the present random number generator.
570
571 Normally, you will create a local generator by one of the constructor calls. But sometimes
572 it is useful to have all program parts access the same generator.
573 */
575 {
576 return global_;
577 }
578
579 static UInt32 factorForUniformInt(UInt32 range)
580 {
581 return (range > 2147483648U || range == 0)
582 ? 1
583 : 2*(2147483648U / ceilPower2(range));
584 }
585
586 static RandomNumberGenerator global_;
587};
588
589template <class Engine>
590RandomNumberGenerator<Engine> RandomNumberGenerator<Engine>::global_(RandomSeed);
591
592
593template <class Engine>
595{
596 if(normalCachedValid_)
597 {
598 normalCachedValid_ = false;
599 return normalCached_;
600 }
601 else
602 {
603 double x1, x2, w;
604 do
605 {
606 x1 = uniform(-1.0, 1.0);
607 x2 = uniform(-1.0, 1.0);
608 w = x1 * x1 + x2 * x2;
609 }
610 while ( w > 1.0 || w == 0.0);
611
612 w = std::sqrt( -2.0 * std::log( w ) / w );
613
614 normalCached_ = x2 * w;
615 normalCachedValid_ = true;
616
617 return x1 * w;
618 }
619}
620
621 /** Shorthand for the TT800 random number generator class.
622 */
624
625 /** Shorthand for the TT800 random number generator class (same as RandomTT800).
626 */
628
629 /** Shorthand for the MT19937 random number generator class.
630 */
632
633 /** Shorthand for the MT19937 random number generator class (same as RandomMT19937).
634 */
636
637 /** Access the global (program-wide) instance of the TT800 random number generator.
638 */
640
641 /** Access the global (program-wide) instance of the MT19937 random number generator.
642 */
644
645template <class Engine>
646class FunctorTraits<RandomNumberGenerator<Engine> >
647{
648 public:
649 typedef RandomNumberGenerator<Engine> type;
650
651 typedef VigraTrueType isInitializer;
652
653 typedef VigraFalseType isUnaryFunctor;
654 typedef VigraFalseType isBinaryFunctor;
655 typedef VigraFalseType isTernaryFunctor;
656
657 typedef VigraFalseType isUnaryAnalyser;
658 typedef VigraFalseType isBinaryAnalyser;
659 typedef VigraFalseType isTernaryAnalyser;
660};
661
662
663/** Functor to create uniformly distributed integer random numbers.
664
665 This functor encapsulates the appropriate functions of the given random number
666 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
667 in an STL-compatible interface.
668
669 <b>Traits defined:</b>
670
671 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim
672 and
673 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor \endverbatim
674 are true (<tt>VigraTrueType</tt>).
675*/
676template <class Engine = MersenneTwister>
678{
679 UInt32 lower_, difference_, factor_;
680 Engine const & generator_;
681 bool useLowBits_;
682
683 public:
684
685 typedef UInt32 argument_type; ///< STL required functor argument type
686 typedef UInt32 result_type; ///< STL required functor result type
687
688 /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>)
689 using the given engine.
690
691 That is, the generated numbers satisfy 0 &lt;= i &lt; 2<sup>32</sup>.
692 */
693 explicit UniformIntRandomFunctor(Engine const & generator = Engine::global() )
694 : lower_(0), difference_(0xffffffff), factor_(1),
695 generator_(generator),
696 useLowBits_(true)
697 {}
698
699 /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>]
700 using the given engine.
701
702 That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
703 \a useLowBits should be set to <tt>false</tt> when the engine generates
704 random numbers whose low bits are significantly less random than the high
705 bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>,
706 but is necessary for simpler linear congruential generators.
707 */
709 Engine const & generator = Engine::global(),
710 bool useLowBits = true)
711 : lower_(lower), difference_(upper-lower),
712 factor_(Engine::factorForUniformInt(difference_ + 1)),
713 generator_(generator),
714 useLowBits_(useLowBits)
715 {
716 vigra_precondition(lower < upper,
717 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
718 }
719
720 /** Return a random number as specified in the constructor.
721 */
723 {
724 if(difference_ == 0xffffffff) // lower_ is necessarily 0
725 return generator_();
726 else if(useLowBits_)
727 return generator_.uniformInt(difference_+1) + lower_;
728 else
729 {
730 UInt32 res = generator_() / factor_;
731
732 // Use rejection method to avoid quantization bias.
733 // On average, we will need two raw random numbers to generate one.
734 while(res > difference_)
735 res = generator_() / factor_;
736 return res + lower_;
737 }
738 }
739
740 /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>).
741
742 That is, 0 &lt;= i &lt; <tt>beyond</tt>. This is a required interface for
743 <tt>std::random_shuffle</tt>. It ignores the limits specified
744 in the constructor and the flag <tt>useLowBits</tt>.
745 */
747 {
748 if(beyond < 2)
749 return 0;
750
751 return generator_.uniformInt(beyond);
752 }
753};
754
755template <class Engine>
756class FunctorTraits<UniformIntRandomFunctor<Engine> >
757{
758 public:
759 typedef UniformIntRandomFunctor<Engine> type;
760
761 typedef VigraTrueType isInitializer;
762
763 typedef VigraTrueType isUnaryFunctor;
764 typedef VigraFalseType isBinaryFunctor;
765 typedef VigraFalseType isTernaryFunctor;
766
767 typedef VigraFalseType isUnaryAnalyser;
768 typedef VigraFalseType isBinaryAnalyser;
769 typedef VigraFalseType isTernaryAnalyser;
770};
771
772/** Functor to create uniformly distributed double-precision random numbers.
773
774 This functor encapsulates the function <tt>uniform()</tt> of the given random number
775 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
776 in an STL-compatible interface.
777
778 <b>Traits defined:</b>
779
780 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim
781 is true (<tt>VigraTrueType</tt>).
782*/
783template <class Engine = MersenneTwister>
785{
786 double offset_, scale_;
787 Engine const & generator_;
788
789 public:
790
791 typedef double result_type; ///< STL required functor result type
792
793 /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0]
794 using the given engine.
795
796 That is, the generated numbers satisfy 0.0 &lt;= i &lt;= 1.0.
797 */
798 UniformRandomFunctor(Engine const & generator = Engine::global())
799 : offset_(0.0),
800 scale_(1.0),
801 generator_(generator)
802 {}
803
804 /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>]
805 using the given engine.
806
807 That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
808 */
809 UniformRandomFunctor(double lower, double upper,
810 Engine & generator = Engine::global())
811 : offset_(lower),
812 scale_(upper - lower),
813 generator_(generator)
814 {
815 vigra_precondition(lower < upper,
816 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
817 }
818
819 /** Return a random number as specified in the constructor.
820 */
821 double operator()() const
822 {
823 return generator_.uniform() * scale_ + offset_;
824 }
825};
826
827template <class Engine>
828class FunctorTraits<UniformRandomFunctor<Engine> >
829{
830 public:
831 typedef UniformRandomFunctor<Engine> type;
832
833 typedef VigraTrueType isInitializer;
834
835 typedef VigraFalseType isUnaryFunctor;
836 typedef VigraFalseType isBinaryFunctor;
837 typedef VigraFalseType isTernaryFunctor;
838
839 typedef VigraFalseType isUnaryAnalyser;
840 typedef VigraFalseType isBinaryAnalyser;
841 typedef VigraFalseType isTernaryAnalyser;
842};
843
844/** Functor to create normal variate random numbers.
845
846 This functor encapsulates the function <tt>normal()</tt> of the given random number
847 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
848 in an STL-compatible interface.
849
850 <b>Traits defined:</b>
851
852 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim
853 is true (<tt>VigraTrueType</tt>).
854*/
855template <class Engine = MersenneTwister>
857{
858 double mean_, stddev_;
859 Engine const & generator_;
860
861 public:
862
863 typedef double result_type; ///< STL required functor result type
864
865 /** Create functor for standard normal random numbers
866 using the given engine.
867
868 That is, mean is 0.0 and standard deviation is 1.0.
869 */
870 NormalRandomFunctor(Engine const & generator = Engine::global())
871 : mean_(0.0),
872 stddev_(1.0),
873 generator_(generator)
874 {}
875
876 /** Create functor for normal random numbers with given mean and standard deviation
877 using the given engine.
878 */
879 NormalRandomFunctor(double mean, double stddev,
880 Engine & generator = Engine::global())
881 : mean_(mean),
882 stddev_(stddev),
883 generator_(generator)
884 {
885 vigra_precondition(stddev > 0.0,
886 "NormalRandomFunctor(): standard deviation must be positive.");
887 }
888
889 /** Return a random number as specified in the constructor.
890 */
891 double operator()() const
892 {
893 return generator_.normal() * stddev_ + mean_;
894 }
895
896};
897
898template <class Engine>
899class FunctorTraits<NormalRandomFunctor<Engine> >
900{
901 public:
902 typedef UniformRandomFunctor<Engine> type;
903
904 typedef VigraTrueType isInitializer;
905
906 typedef VigraFalseType isUnaryFunctor;
907 typedef VigraFalseType isBinaryFunctor;
908 typedef VigraFalseType isTernaryFunctor;
909
910 typedef VigraFalseType isUnaryAnalyser;
911 typedef VigraFalseType isBinaryAnalyser;
912 typedef VigraFalseType isTernaryAnalyser;
913};
914
915//@}
916
917} // namespace vigra
918
919#endif // VIGRA_RANDOM_HXX
NormalRandomFunctor(Engine const &generator=Engine::global())
Definition random.hxx:870
double result_type
STL required functor result type.
Definition random.hxx:863
NormalRandomFunctor(double mean, double stddev, Engine &generator=Engine::global())
Definition random.hxx:879
double operator()() const
Definition random.hxx:891
Definition random.hxx:346
RandomNumberGenerator(Iterator init, UInt32 length)
Definition random.hxx:404
RandomNumberGenerator()
Definition random.hxx:357
double uniform(double lower, double upper) const
Definition random.hxx:544
void seed(Iterator init, UInt32 length)
Definition random.hxx:453
UInt32 operator()() const
Definition random.hxx:463
void seed(RandomSeedTag)
Definition random.hxx:424
UInt32 uniformInt() const
Definition random.hxx:472
UInt32 uniformInt(UInt32 beyond) const
Definition random.hxx:500
void seed(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:438
RandomNumberGenerator(RandomSeedTag)
Definition random.hxx:372
double normal() const
Definition random.hxx:594
double uniform53() const
Definition random.hxx:523
double normal(double mean, double stddev) const
Definition random.hxx:562
RandomNumberGenerator(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:388
static RandomNumberGenerator & global()
Definition random.hxx:574
double uniform() const
Definition random.hxx:534
UniformIntRandomFunctor(UInt32 lower, UInt32 upper, Engine const &generator=Engine::global(), bool useLowBits=true)
Definition random.hxx:708
UInt32 operator()(UInt32 beyond) const
Definition random.hxx:746
UInt32 operator()() const
Definition random.hxx:722
UInt32 argument_type
STL required functor argument type.
Definition random.hxx:685
UniformIntRandomFunctor(Engine const &generator=Engine::global())
Definition random.hxx:693
UInt32 result_type
STL required functor result type.
Definition random.hxx:686
UniformRandomFunctor(double lower, double upper, Engine &generator=Engine::global())
Definition random.hxx:809
double result_type
STL required functor result type.
Definition random.hxx:791
UniformRandomFunctor(Engine const &generator=Engine::global())
Definition random.hxx:798
double operator()() const
Definition random.hxx:821
LookupTag< TAG, A >::result_type get(A const &a)
Definition accumulator.hxx:2942
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition sized_int.hxx:183
RandomNumberGenerator< detail::RandomState< detail::TT800 > > TemperedTwister
Definition random.hxx:627
RandomNumberGenerator< detail::RandomState< detail::MT19937 > > MersenneTwister
Definition random.hxx:635
UInt32 ceilPower2(UInt32 x)
Round up to the nearest power of 2.
Definition mathutil.hxx:294
RandomNumberGenerator< detail::RandomState< detail::MT19937 > > RandomMT19937
Definition random.hxx:631
RandomTT800 & randomTT800()
Definition random.hxx:639
RandomNumberGenerator< detail::RandomState< detail::TT800 > > RandomTT800
Definition random.hxx:623
RandomMT19937 & randomMT19937()
Definition random.hxx:643

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.12.2 (Mon Apr 14 2025)