[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/random.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2008 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 00037 #ifndef VIGRA_RANDOM_HXX 00038 #define VIGRA_RANDOM_HXX 00039 00040 #include "mathutil.hxx" 00041 #include "functortraits.hxx" 00042 #include "array_vector.hxx" 00043 00044 #include <ctime> 00045 00046 // includes to get the current process and thread IDs 00047 // to be used for automated seeding 00048 #ifdef _MSC_VER 00049 #include <vigra/windows.h> // for GetCurrentProcessId() and GetCurrentThreadId() 00050 #endif 00051 00052 #if __linux__ 00053 #include <unistd.h> // for getpid() 00054 #include <sys/syscall.h> // for SYS_gettid 00055 #endif 00056 00057 #if __APPLE__ 00058 #include <unistd.h> // for getpid() 00059 #include <sys/syscall.h> // SYS_thread_selfid 00060 #include <AvailabilityMacros.h> // to check if we are on MacOS 10.6 or later 00061 #endif 00062 00063 namespace vigra { 00064 00065 enum RandomSeedTag { RandomSeed }; 00066 00067 namespace detail { 00068 00069 enum RandomEngineTag { TT800, MT19937 }; 00070 00071 00072 template<RandomEngineTag EngineTag> 00073 struct RandomState; 00074 00075 template <RandomEngineTag EngineTag> 00076 void seed(UInt32 theSeed, RandomState<EngineTag> & engine) 00077 { 00078 engine.state_[0] = theSeed; 00079 for(UInt32 i=1; i<RandomState<EngineTag>::N; ++i) 00080 { 00081 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i; 00082 } 00083 } 00084 00085 template <class Iterator, RandomEngineTag EngineTag> 00086 void seed(Iterator init, UInt32 key_length, RandomState<EngineTag> & engine) 00087 { 00088 const UInt32 N = RandomState<EngineTag>::N; 00089 int k = (int)std::max(N, key_length); 00090 UInt32 i = 1, j = 0; 00091 Iterator data = init; 00092 for (; k; --k) 00093 { 00094 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U)) 00095 + *data + j; /* non linear */ 00096 ++i; ++j; ++data; 00097 00098 if (i >= N) 00099 { 00100 engine.state_[0] = engine.state_[N-1]; 00101 i=1; 00102 } 00103 if (j>=key_length) 00104 { 00105 j=0; 00106 data = init; 00107 } 00108 } 00109 00110 for (k=N-1; k; --k) 00111 { 00112 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U)) 00113 - i; /* non linear */ 00114 ++i; 00115 if (i>=N) 00116 { 00117 engine.state_[0] = engine.state_[N-1]; 00118 i=1; 00119 } 00120 } 00121 00122 engine.state_[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */ 00123 } 00124 00125 template <RandomEngineTag EngineTag> 00126 void seed(RandomSeedTag, RandomState<EngineTag> & engine) 00127 { 00128 static UInt32 globalCount = 0; 00129 ArrayVector<UInt32> seedData; 00130 00131 seedData.push_back((UInt32)time(0)); 00132 seedData.push_back((UInt32)clock()); 00133 seedData.push_back(++globalCount); 00134 00135 std::size_t ptr((char*)&engine - (char*)0); 00136 seedData.push_back((UInt32)(ptr & 0xffffffff)); 00137 seedData.push_back((UInt32)(ptr >> 32)); 00138 00139 #ifdef _MSC_VER 00140 seedData.push_back((UInt32)GetCurrentProcessId()); 00141 seedData.push_back((UInt32)GetCurrentThreadId()); 00142 #endif 00143 00144 #ifdef __linux__ 00145 seedData.push_back((UInt32)getpid()); 00146 # ifdef SYS_gettid 00147 seedData.push_back((UInt32)syscall(SYS_gettid)); 00148 # endif 00149 #endif 00150 00151 #ifdef __APPLE__ 00152 seedData.push_back((UInt32)getpid()); 00153 #if defined(SYS_thread_selfid) && (MAC_OS_X_VERSION_MIN_REQUIRED >= MAC_OS_X_VERSION_10_6) 00154 // SYS_thread_selfid was introduced in MacOS 10.6 00155 seedData.push_back((UInt32)syscall(SYS_thread_selfid)); 00156 #endif 00157 #endif 00158 00159 seed(seedData.begin(), seedData.size(), engine); 00160 } 00161 00162 /* Tempered twister TT800 by M. Matsumoto */ 00163 template<> 00164 struct RandomState<TT800> 00165 { 00166 static const UInt32 N = 25, M = 7; 00167 00168 mutable UInt32 state_[N]; 00169 mutable UInt32 current_; 00170 00171 RandomState() 00172 : current_(0) 00173 { 00174 UInt32 seeds[N] = { 00175 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23, 00176 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825, 00177 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f, 00178 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9, 00179 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb 00180 }; 00181 00182 for(UInt32 i=0; i<N; ++i) 00183 state_[i] = seeds[i]; 00184 } 00185 00186 protected: 00187 00188 UInt32 get() const 00189 { 00190 if(current_ == N) 00191 generateNumbers<void>(); 00192 00193 UInt32 y = state_[current_++]; 00194 y ^= (y << 7) & 0x2b5b2500; 00195 y ^= (y << 15) & 0xdb8b0000; 00196 return y ^ (y >> 16); 00197 } 00198 00199 template <class DUMMY> 00200 void generateNumbers() const; 00201 00202 void seedImpl(RandomSeedTag) 00203 { 00204 seed(RandomSeed, *this); 00205 } 00206 00207 void seedImpl(UInt32 theSeed) 00208 { 00209 seed(theSeed, *this); 00210 } 00211 00212 template<class Iterator> 00213 void seedImpl(Iterator init, UInt32 length) 00214 { 00215 seed(init, length, *this); 00216 } 00217 }; 00218 00219 template <class DUMMY> 00220 void RandomState<TT800>::generateNumbers() const 00221 { 00222 UInt32 mag01[2]= { 0x0, 0x8ebfd028 }; 00223 00224 for(UInt32 i=0; i<N-M; ++i) 00225 { 00226 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2]; 00227 } 00228 for (UInt32 i=N-M; i<N; ++i) 00229 { 00230 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2]; 00231 } 00232 current_ = 0; 00233 } 00234 00235 /* Mersenne twister MT19937 by M. Matsumoto */ 00236 template<> 00237 struct RandomState<MT19937> 00238 { 00239 static const UInt32 N = 624, M = 397; 00240 00241 mutable UInt32 state_[N]; 00242 mutable UInt32 current_; 00243 00244 RandomState() 00245 : current_(0) 00246 { 00247 seed(19650218U, *this); 00248 } 00249 00250 protected: 00251 00252 UInt32 get() const 00253 { 00254 if(current_ == N) 00255 generateNumbers<void>(); 00256 00257 UInt32 x = state_[current_++]; 00258 x ^= (x >> 11); 00259 x ^= (x << 7) & 0x9D2C5680U; 00260 x ^= (x << 15) & 0xEFC60000U; 00261 return x ^ (x >> 18); 00262 } 00263 00264 template <class DUMMY> 00265 void generateNumbers() const; 00266 00267 static UInt32 twiddle(UInt32 u, UInt32 v) 00268 { 00269 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1) 00270 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U); 00271 } 00272 00273 void seedImpl(RandomSeedTag) 00274 { 00275 seed(RandomSeed, *this); 00276 generateNumbers<void>(); 00277 } 00278 00279 void seedImpl(UInt32 theSeed) 00280 { 00281 seed(theSeed, *this); 00282 generateNumbers<void>(); 00283 } 00284 00285 template<class Iterator> 00286 void seedImpl(Iterator init, UInt32 length) 00287 { 00288 seed(19650218U, *this); 00289 seed(init, length, *this); 00290 generateNumbers<void>(); 00291 } 00292 }; 00293 00294 template <class DUMMY> 00295 void RandomState<MT19937>::generateNumbers() const 00296 { 00297 for (unsigned int i = 0; i < (N - M); ++i) 00298 { 00299 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]); 00300 } 00301 for (unsigned int i = N - M; i < (N - 1); ++i) 00302 { 00303 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]); 00304 } 00305 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]); 00306 current_ = 0; 00307 } 00308 00309 } // namespace detail 00310 00311 00312 /** \addtogroup RandomNumberGeneration Random Number Generation 00313 00314 High-quality random number generators and related functors. 00315 */ 00316 //@{ 00317 00318 /** Generic random number generator. 00319 00320 The actual generator is passed in the template argument <tt>Engine</tt>. Two generators 00321 are currently available: 00322 <ul> 00323 <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. 00324 <li> <tt>RandomTT800</tt>: (default) The Tempered Twister, a simpler predecessor of the Mersenne Twister with period length 2<sup>800</sup>. 00325 </ul> 00326 00327 Both generators have been designed by <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html">Makoto Matsumoto</a>. 00328 00329 <b>Traits defined:</b> 00330 00331 \verbatim FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer \endverbatim 00332 is true (<tt>VigraTrueType</tt>). 00333 */ 00334 template <class Engine = detail::RandomState<detail::TT800> > 00335 class RandomNumberGenerator 00336 : public Engine 00337 { 00338 mutable double normalCached_; 00339 mutable bool normalCachedValid_; 00340 00341 public: 00342 00343 /** Create a new random generator object with standard seed. 00344 00345 Due to standard seeding, the random numbers generated will always be the same. 00346 This is useful for debugging. 00347 */ 00348 RandomNumberGenerator() 00349 : normalCached_(0.0), 00350 normalCachedValid_(false) 00351 {} 00352 00353 /** Create a new random generator object with a random seed. 00354 00355 The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt> 00356 values. 00357 00358 <b>Usage:</b> 00359 \code 00360 RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed); 00361 \endcode 00362 */ 00363 RandomNumberGenerator(RandomSeedTag) 00364 : normalCached_(0.0), 00365 normalCachedValid_(false) 00366 { 00367 this->seedImpl(RandomSeed); 00368 } 00369 00370 /** Create a new random generator object from the given seed. 00371 00372 The same seed will always produce identical random sequences. 00373 */ 00374 RandomNumberGenerator(UInt32 theSeed) 00375 : normalCached_(0.0), 00376 normalCachedValid_(false) 00377 { 00378 this->seedImpl(theSeed); 00379 } 00380 00381 /** Create a new random generator object from the given seed sequence. 00382 00383 Longer seed sequences lead to better initialization in the sense that the generator's 00384 state space is covered much better than is possible with 32-bit seeds alone. 00385 */ 00386 template<class Iterator> 00387 RandomNumberGenerator(Iterator init, UInt32 length) 00388 : normalCached_(0.0), 00389 normalCachedValid_(false) 00390 { 00391 this->seedImpl(init, length); 00392 } 00393 00394 00395 /** Re-initialize the random generator object with a random seed. 00396 00397 The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt> 00398 values. 00399 00400 <b>Usage:</b> 00401 \code 00402 RandomNumberGenerator<> rnd = ...; 00403 ... 00404 rnd.seed(RandomSeed); 00405 \endcode 00406 */ 00407 void seed(RandomSeedTag) 00408 { 00409 this->seedImpl(RandomSeed); 00410 normalCachedValid_ = false; 00411 } 00412 00413 /** Re-initialize the random generator object from the given seed. 00414 00415 The same seed will always produce identical random sequences. 00416 */ 00417 void seed(UInt32 theSeed) 00418 { 00419 this->seedImpl(theSeed); 00420 normalCachedValid_ = false; 00421 } 00422 00423 /** Re-initialize the random generator object from the given seed sequence. 00424 00425 Longer seed sequences lead to better initialization in the sense that the generator's 00426 state space is covered much better than is possible with 32-bit seeds alone. 00427 */ 00428 template<class Iterator> 00429 void seed(Iterator init, UInt32 length) 00430 { 00431 this->seedImpl(init, length); 00432 normalCachedValid_ = false; 00433 } 00434 00435 /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>). 00436 00437 That is, 0 <= i < 2<sup>32</sup>. 00438 */ 00439 UInt32 operator()() const 00440 { 00441 return this->get(); 00442 } 00443 00444 /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>). 00445 00446 That is, 0 <= i < 2<sup>32</sup>. 00447 */ 00448 UInt32 uniformInt() const 00449 { 00450 return this->get(); 00451 } 00452 00453 00454 #if 0 // difficult implementation necessary if low bits are not sufficiently random 00455 // in [0,beyond) 00456 UInt32 uniformInt(UInt32 beyond) const 00457 { 00458 if(beyond < 2) 00459 return 0; 00460 00461 UInt32 factor = factorForUniformInt(beyond); 00462 UInt32 res = this->get() / factor; 00463 00464 // Use rejection method to avoid quantization bias. 00465 // On average, we will need two raw random numbers to generate one. 00466 while(res >= beyond) 00467 res = this->get() / factor; 00468 return res; 00469 } 00470 #endif /* #if 0 */ 00471 00472 /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>). 00473 00474 That is, 0 <= i < <tt>beyond</tt>. 00475 */ 00476 UInt32 uniformInt(UInt32 beyond) const 00477 { 00478 // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is 00479 // the case for TT800 and MT19937 00480 if(beyond < 2) 00481 return 0; 00482 00483 UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond; 00484 UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder; 00485 UInt32 res = this->get(); 00486 00487 // Use rejection method to avoid quantization bias. 00488 // We will need two raw random numbers in amortized worst case. 00489 while(res > lastSafeValue) 00490 res = this->get(); 00491 return res % beyond; 00492 } 00493 00494 /** Return a uniformly distributed double-precision random number in [0.0, 1.0). 00495 00496 That is, 0.0 <= i < 1.0. All 53-bit bits of the mantissa are random (two 32-bit integers are used to 00497 create this number). 00498 */ 00499 double uniform53() const 00500 { 00501 // make full use of the entire 53-bit mantissa of a double, by Isaku Wada 00502 return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0); 00503 } 00504 00505 /** Return a uniformly distributed double-precision random number in [0.0, 1.0]. 00506 00507 That is, 0.0 <= i <= 1.0. This number is computed by <tt>uniformInt()</tt> / (2<sup>32</sup> - 1), 00508 so it has effectively only 32 random bits. 00509 */ 00510 double uniform() const 00511 { 00512 return (double)this->get() / 4294967295.0; 00513 } 00514 00515 /** Return a uniformly distributed double-precision random number in [lower, upper]. 00516 00517 That is, <tt>lower</tt> <= i <= <tt>upper</tt>. This number is computed 00518 from <tt>uniform()</tt>, so it has effectively only 32 random bits. 00519 */ 00520 double uniform(double lower, double upper) const 00521 { 00522 vigra_precondition(lower < upper, 00523 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound."); 00524 return uniform() * (upper-lower) + lower; 00525 } 00526 00527 /** Return a standard normal variate (Gaussian) random number. 00528 00529 Mean is zero, standard deviation is 1.0. It uses the polar form of the 00530 Box-Muller transform. 00531 */ 00532 double normal() const; 00533 00534 /** Return a normal variate (Gaussian) random number with the given mean and standard deviation. 00535 00536 It uses the polar form of the Box-Muller transform. 00537 */ 00538 double normal(double mean, double stddev) const 00539 { 00540 vigra_precondition(stddev > 0.0, 00541 "RandomNumberGenerator::normal(): standard deviation must be positive."); 00542 return normal()*stddev + mean; 00543 } 00544 00545 /** Access the global (program-wide) instance of the present random number generator. 00546 00547 Normally, you will create a local generator by one of the constructor calls. But sometimes 00548 it is useful to have all program parts access the same generator. 00549 */ 00550 static RandomNumberGenerator & global() 00551 { 00552 static RandomNumberGenerator generator; 00553 return generator; 00554 } 00555 00556 static UInt32 factorForUniformInt(UInt32 range) 00557 { 00558 return (range > 2147483648U || range == 0) 00559 ? 1 00560 : 2*(2147483648U / ceilPower2(range)); 00561 } 00562 }; 00563 00564 template <class Engine> 00565 double RandomNumberGenerator<Engine>::normal() const 00566 { 00567 if(normalCachedValid_) 00568 { 00569 normalCachedValid_ = false; 00570 return normalCached_; 00571 } 00572 else 00573 { 00574 double x1, x2, w; 00575 do 00576 { 00577 x1 = uniform(-1.0, 1.0); 00578 x2 = uniform(-1.0, 1.0); 00579 w = x1 * x1 + x2 * x2; 00580 } 00581 while ( w > 1.0 || w == 0.0); 00582 00583 w = std::sqrt( -2.0 * std::log( w ) / w ); 00584 00585 normalCached_ = x2 * w; 00586 normalCachedValid_ = true; 00587 00588 return x1 * w; 00589 } 00590 } 00591 00592 /** Shorthand for the TT800 random number generator class. 00593 */ 00594 typedef RandomNumberGenerator<> RandomTT800; 00595 00596 /** Shorthand for the TT800 random number generator class (same as RandomTT800). 00597 */ 00598 typedef RandomNumberGenerator<> TemperedTwister; 00599 00600 /** Shorthand for the MT19937 random number generator class. 00601 */ 00602 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > RandomMT19937; 00603 00604 /** Shorthand for the MT19937 random number generator class (same as RandomMT19937). 00605 */ 00606 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > MersenneTwister; 00607 00608 /** Access the global (program-wide) instance of the TT800 random number generator. 00609 */ 00610 inline RandomTT800 & randomTT800() { return RandomTT800::global(); } 00611 00612 /** Access the global (program-wide) instance of the MT19937 random number generator. 00613 */ 00614 inline RandomMT19937 & randomMT19937() { return RandomMT19937::global(); } 00615 00616 template <class Engine> 00617 class FunctorTraits<RandomNumberGenerator<Engine> > 00618 { 00619 public: 00620 typedef RandomNumberGenerator<Engine> type; 00621 00622 typedef VigraTrueType isInitializer; 00623 00624 typedef VigraFalseType isUnaryFunctor; 00625 typedef VigraFalseType isBinaryFunctor; 00626 typedef VigraFalseType isTernaryFunctor; 00627 00628 typedef VigraFalseType isUnaryAnalyser; 00629 typedef VigraFalseType isBinaryAnalyser; 00630 typedef VigraFalseType isTernaryAnalyser; 00631 }; 00632 00633 00634 /** Functor to create uniformly distributed integer random numbers. 00635 00636 This functor encapsulates the appropriate functions of the given random number 00637 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>) 00638 in an STL-compatible interface. 00639 00640 <b>Traits defined:</b> 00641 00642 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim 00643 and 00644 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor \endverbatim 00645 are true (<tt>VigraTrueType</tt>). 00646 */ 00647 template <class Engine = RandomTT800> 00648 class UniformIntRandomFunctor 00649 { 00650 UInt32 lower_, difference_, factor_; 00651 Engine const & generator_; 00652 bool useLowBits_; 00653 00654 public: 00655 00656 typedef UInt32 argument_type; ///< STL required functor argument type 00657 typedef UInt32 result_type; ///< STL required functor result type 00658 00659 /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>) 00660 using the given engine. 00661 00662 That is, the generated numbers satisfy 0 <= i < 2<sup>32</sup>. 00663 */ 00664 explicit UniformIntRandomFunctor(Engine const & generator = Engine::global() ) 00665 : lower_(0), difference_(0xffffffff), factor_(1), 00666 generator_(generator), 00667 useLowBits_(true) 00668 {} 00669 00670 /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>] 00671 using the given engine. 00672 00673 That is, the generated numbers satisfy <tt>lower</tt> <= i <= <tt>upper</tt>. 00674 \a useLowBits should be set to <tt>false</tt> when the engine generates 00675 random numbers whose low bits are significantly less random than the high 00676 bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>, 00677 but is necessary for simpler linear congruential generators. 00678 */ 00679 UniformIntRandomFunctor(UInt32 lower, UInt32 upper, 00680 Engine const & generator = Engine::global(), 00681 bool useLowBits = true) 00682 : lower_(lower), difference_(upper-lower), 00683 factor_(Engine::factorForUniformInt(difference_ + 1)), 00684 generator_(generator), 00685 useLowBits_(useLowBits) 00686 { 00687 vigra_precondition(lower < upper, 00688 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound."); 00689 } 00690 00691 /** Return a random number as specified in the constructor. 00692 */ 00693 UInt32 operator()() const 00694 { 00695 if(difference_ == 0xffffffff) // lower_ is necessarily 0 00696 return generator_(); 00697 else if(useLowBits_) 00698 return generator_.uniformInt(difference_+1) + lower_; 00699 else 00700 { 00701 UInt32 res = generator_() / factor_; 00702 00703 // Use rejection method to avoid quantization bias. 00704 // On average, we will need two raw random numbers to generate one. 00705 while(res > difference_) 00706 res = generator_() / factor_; 00707 return res + lower_; 00708 } 00709 } 00710 00711 /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>). 00712 00713 That is, 0 <= i < <tt>beyond</tt>. This is a required interface for 00714 <tt>std::random_shuffle</tt>. It ignores the limits specified 00715 in the constructor and the flag <tt>useLowBits</tt>. 00716 */ 00717 UInt32 operator()(UInt32 beyond) const 00718 { 00719 if(beyond < 2) 00720 return 0; 00721 00722 return generator_.uniformInt(beyond); 00723 } 00724 }; 00725 00726 template <class Engine> 00727 class FunctorTraits<UniformIntRandomFunctor<Engine> > 00728 { 00729 public: 00730 typedef UniformIntRandomFunctor<Engine> type; 00731 00732 typedef VigraTrueType isInitializer; 00733 00734 typedef VigraTrueType isUnaryFunctor; 00735 typedef VigraFalseType isBinaryFunctor; 00736 typedef VigraFalseType isTernaryFunctor; 00737 00738 typedef VigraFalseType isUnaryAnalyser; 00739 typedef VigraFalseType isBinaryAnalyser; 00740 typedef VigraFalseType isTernaryAnalyser; 00741 }; 00742 00743 /** Functor to create uniformly distributed double-precision random numbers. 00744 00745 This functor encapsulates the function <tt>uniform()</tt> of the given random number 00746 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>) 00747 in an STL-compatible interface. 00748 00749 <b>Traits defined:</b> 00750 00751 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim 00752 is true (<tt>VigraTrueType</tt>). 00753 */ 00754 template <class Engine = RandomTT800> 00755 class UniformRandomFunctor 00756 { 00757 double offset_, scale_; 00758 Engine const & generator_; 00759 00760 public: 00761 00762 typedef double result_type; ///< STL required functor result type 00763 00764 /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0] 00765 using the given engine. 00766 00767 That is, the generated numbers satisfy 0.0 <= i <= 1.0. 00768 */ 00769 UniformRandomFunctor(Engine const & generator = Engine::global()) 00770 : offset_(0.0), 00771 scale_(1.0), 00772 generator_(generator) 00773 {} 00774 00775 /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>] 00776 using the given engine. 00777 00778 That is, the generated numbers satisfy <tt>lower</tt> <= i <= <tt>upper</tt>. 00779 */ 00780 UniformRandomFunctor(double lower, double upper, 00781 Engine & generator = Engine::global()) 00782 : offset_(lower), 00783 scale_(upper - lower), 00784 generator_(generator) 00785 { 00786 vigra_precondition(lower < upper, 00787 "UniformRandomFunctor(): lower bound must be smaller than upper bound."); 00788 } 00789 00790 /** Return a random number as specified in the constructor. 00791 */ 00792 double operator()() const 00793 { 00794 return generator_.uniform() * scale_ + offset_; 00795 } 00796 }; 00797 00798 template <class Engine> 00799 class FunctorTraits<UniformRandomFunctor<Engine> > 00800 { 00801 public: 00802 typedef UniformRandomFunctor<Engine> type; 00803 00804 typedef VigraTrueType isInitializer; 00805 00806 typedef VigraFalseType isUnaryFunctor; 00807 typedef VigraFalseType isBinaryFunctor; 00808 typedef VigraFalseType isTernaryFunctor; 00809 00810 typedef VigraFalseType isUnaryAnalyser; 00811 typedef VigraFalseType isBinaryAnalyser; 00812 typedef VigraFalseType isTernaryAnalyser; 00813 }; 00814 00815 /** Functor to create normal variate random numbers. 00816 00817 This functor encapsulates the function <tt>normal()</tt> of the given random number 00818 <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>) 00819 in an STL-compatible interface. 00820 00821 <b>Traits defined:</b> 00822 00823 \verbatim FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer \endverbatim 00824 is true (<tt>VigraTrueType</tt>). 00825 */ 00826 template <class Engine = RandomTT800> 00827 class NormalRandomFunctor 00828 { 00829 double mean_, stddev_; 00830 Engine const & generator_; 00831 00832 public: 00833 00834 typedef double result_type; ///< STL required functor result type 00835 00836 /** Create functor for standard normal random numbers 00837 using the given engine. 00838 00839 That is, mean is 0.0 and standard deviation is 1.0. 00840 */ 00841 NormalRandomFunctor(Engine const & generator = Engine::global()) 00842 : mean_(0.0), 00843 stddev_(1.0), 00844 generator_(generator) 00845 {} 00846 00847 /** Create functor for normal random numbers with given mean and standard deviation 00848 using the given engine. 00849 */ 00850 NormalRandomFunctor(double mean, double stddev, 00851 Engine & generator = Engine::global()) 00852 : mean_(mean), 00853 stddev_(stddev), 00854 generator_(generator) 00855 { 00856 vigra_precondition(stddev > 0.0, 00857 "NormalRandomFunctor(): standard deviation must be positive."); 00858 } 00859 00860 /** Return a random number as specified in the constructor. 00861 */ 00862 double operator()() const 00863 { 00864 return generator_.normal() * stddev_ + mean_; 00865 } 00866 00867 }; 00868 00869 template <class Engine> 00870 class FunctorTraits<NormalRandomFunctor<Engine> > 00871 { 00872 public: 00873 typedef UniformRandomFunctor<Engine> type; 00874 00875 typedef VigraTrueType isInitializer; 00876 00877 typedef VigraFalseType isUnaryFunctor; 00878 typedef VigraFalseType isBinaryFunctor; 00879 typedef VigraFalseType isTernaryFunctor; 00880 00881 typedef VigraFalseType isUnaryAnalyser; 00882 typedef VigraFalseType isBinaryAnalyser; 00883 typedef VigraFalseType isTernaryAnalyser; 00884 }; 00885 00886 //@} 00887 00888 } // namespace vigra 00889 00890 #endif // VIGRA_RANDOM_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|