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