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

vigra/random.hxx VIGRA

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 &lt;= i &lt; 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 &lt;= i &lt; 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 &lt;= i &lt; <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 &lt;= i &lt; 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 &lt;= i &lt;= 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> &lt;= i &lt;= <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 &lt;= i &lt; 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> &lt;= i &lt;= <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 &lt;= i &lt; <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 &lt;= i &lt;= 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> &lt;= i &lt;= <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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)