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