[ 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     <tt>FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00285 */
00286 template <class Engine = detail::RandomState<detail::TT800> >
00287 class RandomNumberGenerator
00288 : public Engine
00289 {
00290     mutable double normalCached_;
00291     mutable bool normalCachedValid_;
00292     
00293   public:
00294   
00295         /** Create a new random generator object with standard seed.
00296             
00297             Due to standard seeding, the random numbers generated will always be the same. 
00298             This is useful for debugging.
00299         */
00300     RandomNumberGenerator()
00301     : normalCached_(0.0),
00302       normalCachedValid_(false)
00303     {}
00304   
00305         /** Create a new random generator object with a random seed.
00306         
00307             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
00308             values.
00309         
00310             <b>Usage:</b>
00311             \code
00312             RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed);
00313             \endcode
00314         */
00315     RandomNumberGenerator(RandomSeedTag)
00316     : normalCached_(0.0),
00317       normalCachedValid_(false)
00318     {
00319         this->seedImpl(RandomSeed);
00320     }
00321   
00322         /** Create a new random generator object from the given seed.
00323             
00324             The same seed will always produce identical random sequences.
00325         */
00326     RandomNumberGenerator(UInt32 theSeed)
00327     : normalCached_(0.0),
00328       normalCachedValid_(false)
00329     {
00330         this->seedImpl(theSeed);
00331     }
00332 
00333         /** Create a new random generator object from the given seed sequence.
00334             
00335             Longer seed sequences lead to better initialization in the sense that the generator's 
00336             state space is covered much better than is possible with 32-bit seeds alone.
00337         */
00338     template<class Iterator>
00339     RandomNumberGenerator(Iterator init, UInt32 length)
00340     : normalCached_(0.0),
00341       normalCachedValid_(false)
00342     {
00343         this->seedImpl(init, length);
00344     }
00345 
00346   
00347         /** Re-initialize the random generator object with a random seed.
00348         
00349             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
00350             values.
00351         
00352             <b>Usage:</b>
00353             \code
00354             RandomNumberGenerator<> rnd = ...;
00355             ...
00356             rnd.seed(RandomSeed);
00357             \endcode
00358         */
00359     void seed(RandomSeedTag)
00360     {
00361         this->seedImpl(RandomSeed);
00362         normalCachedValid_ = false;
00363     }
00364 
00365         /** Re-initialize the random generator object from the given seed.
00366             
00367             The same seed will always produce identical random sequences.
00368         */
00369     void seed(UInt32 theSeed)
00370     {
00371         this->seedImpl(theSeed);
00372         normalCachedValid_ = false;
00373     }
00374 
00375         /** Re-initialize the random generator object from the given seed sequence.
00376             
00377             Longer seed sequences lead to better initialization in the sense that the generator's 
00378             state space is covered much better than is possible with 32-bit seeds alone.
00379         */
00380     template<class Iterator>
00381     void seed(Iterator init, UInt32 length)
00382     {
00383         this->seedImpl(init, length);
00384         normalCachedValid_ = false;
00385     }
00386 
00387         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
00388             
00389             That is, 0 &lt;= i &lt; 2<sup>32</sup>. 
00390         */
00391     UInt32 operator()() const
00392     {
00393         return this->get();
00394     }
00395 
00396         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
00397             
00398             That is, 0 &lt;= i &lt; 2<sup>32</sup>. 
00399         */
00400     UInt32 uniformInt() const
00401     {
00402         return this->get();
00403     }
00404 
00405 
00406 #if 0 // difficult implementation necessary if low bits are not sufficiently random
00407         // in [0,beyond)
00408     UInt32 uniformInt(UInt32 beyond) const
00409     {
00410         if(beyond < 2)
00411             return 0;
00412 
00413         UInt32 factor = factorForUniformInt(beyond);
00414         UInt32 res = this->get() / factor;
00415 
00416         // Use rejection method to avoid quantization bias.
00417         // On average, we will need two raw random numbers to generate one.
00418         while(res >= beyond)
00419             res = this->get() / factor;
00420         return res;
00421     }
00422 #endif /* #if 0 */
00423 
00424         /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>).
00425             
00426             That is, 0 &lt;= i &lt; <tt>beyond</tt>. 
00427         */
00428     UInt32 uniformInt(UInt32 beyond) const
00429     {
00430         // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is
00431         // the case for TT800 and MT19937
00432         if(beyond < 2)
00433             return 0;
00434 
00435         UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
00436         UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
00437         UInt32 res = this->get();
00438 
00439         // Use rejection method to avoid quantization bias.
00440         // We will need two raw random numbers in amortized worst case.
00441         while(res > lastSafeValue)
00442             res = this->get();
00443         return res % beyond;
00444     }
00445     
00446         /** Return a uniformly distributed double-precision random number in [0.0, 1.0).
00447             
00448             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 
00449             create this number).
00450         */
00451     double uniform53() const
00452     {
00453         // make full use of the entire 53-bit mantissa of a double, by Isaku Wada
00454         return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0); 
00455     }
00456     
00457         /** Return a uniformly distributed double-precision random number in [0.0, 1.0].
00458             
00459             That is, 0.0 &lt;= i &lt;= 1.0. This nuber is computed by <tt>uniformInt()</tt> / 2<sup>32</sup>, 
00460             so it has effectively only 32 random bits. 
00461         */
00462     double uniform() const
00463     {
00464         return (double)this->get() / 4294967295.0;
00465     }
00466 
00467         /** Return a uniformly distributed double-precision random number in [lower, upper].
00468            
00469             That is, <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>. This number is computed 
00470             from <tt>uniform()</tt>, so it has effectively only 32 random bits. 
00471         */
00472     double uniform(double lower, double upper) const
00473     {
00474         vigra_precondition(lower < upper,
00475           "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound."); 
00476         return uniform() * (upper-lower) + lower;
00477     }
00478 
00479         /** Return a standard normal variate (Gaussian) random number.
00480            
00481             Mean is zero, standard deviation is 1.0. It uses the polar form of the 
00482             Box-Muller transform.
00483         */
00484     double normal() const;
00485     
00486         /** Return a normal variate (Gaussian) random number with the given mean and standard deviation.
00487            
00488             It uses the polar form of the Box-Muller transform.
00489         */
00490     double normal(double mean, double stddev) const
00491     {
00492         vigra_precondition(stddev > 0.0,
00493           "RandomNumberGenerator::normal(): standard deviation must be positive."); 
00494         return normal()*stddev + mean;
00495     }
00496     
00497         /** Access the global (program-wide) instance of the present random number generator.
00498         
00499             Normally, you will create a local generator by one of the constructor calls. But sometimes
00500             it is useful to have all program parts access the same generator.
00501         */
00502     static RandomNumberGenerator & global()
00503     {
00504         static RandomNumberGenerator generator;
00505         return generator;
00506     }
00507 
00508     static UInt32 factorForUniformInt(UInt32 range)
00509     {
00510         return (range > 2147483648U || range == 0)
00511                      ? 1
00512                      : 2*(2147483648U / ceilPower2(range));
00513     }
00514 };
00515 
00516 template <class Engine>
00517 double RandomNumberGenerator<Engine>::normal() const
00518 {
00519     if(normalCachedValid_)
00520     {
00521         normalCachedValid_ = false;
00522         return normalCached_;
00523     }
00524     else
00525     {
00526         double x1, x2, w;
00527         do 
00528         {
00529              x1 = uniform(-1.0, 1.0);
00530              x2 = uniform(-1.0, 1.0);
00531              w = x1 * x1 + x2 * x2;
00532         } 
00533         while ( w > 1.0 || w == 0.0);
00534         
00535         w = std::sqrt( -2.0 * std::log( w )  / w );
00536 
00537         normalCached_ = x2 * w;
00538         normalCachedValid_ = true;
00539 
00540         return x1 * w;
00541     }
00542 }
00543 
00544     /** Shorthand for the TT800 random number generator class.
00545     */
00546 typedef RandomNumberGenerator<>  RandomTT800; 
00547 
00548     /** Shorthand for the TT800 random number generator class (same as RandomTT800).
00549     */
00550 typedef RandomNumberGenerator<>  TemperedTwister; 
00551 
00552     /** Shorthand for the MT19937 random number generator class.
00553     */
00554 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > RandomMT19937;
00555 
00556     /** Shorthand for the MT19937 random number generator class (same as RandomMT19937).
00557     */
00558 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > MersenneTwister;
00559 
00560     /** Access the global (program-wide) instance of the TT800 random number generator.
00561     */
00562 inline RandomTT800   & randomTT800()   { return RandomTT800::global(); }
00563 
00564     /** Access the global (program-wide) instance of the MT19937 random number generator.
00565     */
00566 inline RandomMT19937 & randomMT19937() { return RandomMT19937::global(); }
00567 
00568 template <class Engine>
00569 class FunctorTraits<RandomNumberGenerator<Engine> >
00570 {
00571   public:
00572     typedef RandomNumberGenerator<Engine> type;
00573     
00574     typedef VigraTrueType  isInitializer;
00575     
00576     typedef VigraFalseType isUnaryFunctor;
00577     typedef VigraFalseType isBinaryFunctor;
00578     typedef VigraFalseType isTernaryFunctor;
00579     
00580     typedef VigraFalseType isUnaryAnalyser;
00581     typedef VigraFalseType isBinaryAnalyser;
00582     typedef VigraFalseType isTernaryAnalyser;
00583 };
00584 
00585 
00586 /** Functor to create uniformely distributed integer random numbers.
00587 
00588     This functor encapsulates the appropriate functions of the given random number
00589     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00590     in an STL-compatible interface. 
00591     
00592     <b>Traits defined:</b>
00593     
00594     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> and
00595     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor</tt> are true (<tt>VigraTrueType</tt>).
00596 */
00597 template <class Engine = RandomTT800>
00598 class UniformIntRandomFunctor
00599 {
00600     UInt32 lower_, difference_, factor_;
00601     Engine const & generator_;
00602     bool useLowBits_;
00603 
00604   public:
00605   
00606     typedef UInt32 argument_type; ///< STL required functor argument type
00607     typedef UInt32 result_type; ///< STL required functor result type
00608 
00609         /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>) 
00610             using the given engine.
00611             
00612             That is, the generated numbers satisfy 0 &lt;= i &lt; 2<sup>32</sup>.
00613         */
00614     explicit UniformIntRandomFunctor(Engine const & generator = Engine::global() )
00615     : lower_(0), difference_(0xffffffff), factor_(1),
00616       generator_(generator),
00617       useLowBits_(true)
00618     {}
00619     
00620         /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>]
00621             using the given engine.
00622             
00623             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
00624             \a useLowBits should be set to <tt>false</tt> when the engine generates
00625             random numbers whose low bits are significantly less random than the high
00626             bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>,
00627             but is necessary for simpler linear congruential generators.
00628         */
00629     UniformIntRandomFunctor(UInt32 lower, UInt32 upper, 
00630                             Engine const & generator = Engine::global(),
00631                             bool useLowBits = true)
00632     : lower_(lower), difference_(upper-lower), 
00633       factor_(Engine::factorForUniformInt(difference_ + 1)),
00634       generator_(generator),
00635       useLowBits_(useLowBits)
00636     {
00637         vigra_precondition(lower < upper,
00638           "UniformIntRandomFunctor(): lower bound must be smaller than upper bound."); 
00639     }
00640     
00641         /** Return a random number as specified in the constructor.
00642         */
00643     UInt32 operator()() const
00644     {
00645         if(difference_ == 0xffffffff) // lower_ is necessarily 0
00646             return generator_();
00647         else if(useLowBits_)
00648             return generator_.uniformInt(difference_+1) + lower_;
00649         else
00650         {
00651             UInt32 res = generator_() / factor_;
00652 
00653             // Use rejection method to avoid quantization bias.
00654             // On average, we will need two raw random numbers to generate one.
00655             while(res > difference_)
00656                 res = generator_() / factor_;
00657             return res + lower_;
00658         }
00659     }
00660 
00661         /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>).
00662         
00663             That is, 0 &lt;= i &lt; <tt>beyond</tt>. This is a required interface for 
00664             <tt>std::random_shuffle</tt>. It ignores the limits specified 
00665             in the constructor and the flag <tt>useLowBits</tt>.
00666         */
00667     UInt32 operator()(UInt32 beyond) const
00668     {
00669         if(beyond < 2)
00670             return 0;
00671 
00672         return generator_.uniformInt(beyond);
00673     }
00674 };
00675 
00676 template <class Engine>
00677 class FunctorTraits<UniformIntRandomFunctor<Engine> >
00678 {
00679   public:
00680     typedef UniformIntRandomFunctor<Engine> type;
00681     
00682     typedef VigraTrueType  isInitializer;
00683     
00684     typedef VigraTrueType  isUnaryFunctor;
00685     typedef VigraFalseType isBinaryFunctor;
00686     typedef VigraFalseType isTernaryFunctor;
00687     
00688     typedef VigraFalseType isUnaryAnalyser;
00689     typedef VigraFalseType isBinaryAnalyser;
00690     typedef VigraFalseType isTernaryAnalyser;
00691 };
00692 
00693 /** Functor to create uniformely distributed double-precision random numbers.
00694 
00695     This functor encapsulates the function <tt>uniform()</tt> of the given random number
00696     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00697     in an STL-compatible interface. 
00698     
00699     <b>Traits defined:</b>
00700     
00701     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00702 */
00703 template <class Engine = RandomTT800>
00704 class UniformRandomFunctor
00705 {
00706     double offset_, scale_;
00707     Engine const & generator_;
00708 
00709   public:
00710   
00711     typedef double result_type; ///< STL required functor result type
00712 
00713         /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0] 
00714             using the given engine.
00715             
00716             That is, the generated numbers satisfy 0.0 &lt;= i &lt;= 1.0.
00717         */
00718     UniformRandomFunctor(Engine const & generator = Engine::global())
00719     : offset_(0.0),
00720       scale_(1.0),
00721       generator_(generator)
00722     {}
00723 
00724         /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>]
00725             using the given engine.
00726             
00727             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
00728         */
00729     UniformRandomFunctor(double lower, double upper, 
00730                          Engine & generator = Engine::global())
00731     : offset_(lower),
00732       scale_(upper - lower),
00733       generator_(generator)
00734     {
00735         vigra_precondition(lower < upper,
00736           "UniformRandomFunctor(): lower bound must be smaller than upper bound."); 
00737     }
00738     
00739         /** Return a random number as specified in the constructor.
00740         */
00741     double operator()() const
00742     {
00743         return generator_.uniform() * scale_ + offset_;
00744     }
00745 };
00746 
00747 template <class Engine>
00748 class FunctorTraits<UniformRandomFunctor<Engine> >
00749 {
00750   public:
00751     typedef UniformRandomFunctor<Engine> type;
00752     
00753     typedef VigraTrueType  isInitializer;
00754     
00755     typedef VigraFalseType isUnaryFunctor;
00756     typedef VigraFalseType isBinaryFunctor;
00757     typedef VigraFalseType isTernaryFunctor;
00758     
00759     typedef VigraFalseType isUnaryAnalyser;
00760     typedef VigraFalseType isBinaryAnalyser;
00761     typedef VigraFalseType isTernaryAnalyser;
00762 };
00763 
00764 /** Functor to create normal variate random numbers.
00765 
00766     This functor encapsulates the function <tt>normal()</tt> of the given random number
00767     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00768     in an STL-compatible interface. 
00769     
00770     <b>Traits defined:</b>
00771     
00772     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00773 */
00774 template <class Engine = RandomTT800>
00775 class NormalRandomFunctor
00776 {
00777     double mean_, stddev_;
00778     Engine const & generator_;
00779 
00780   public:
00781   
00782     typedef double result_type; ///< STL required functor result type
00783 
00784         /** Create functor for standard normal random numbers 
00785             using the given engine.
00786             
00787             That is, mean is 0.0 and standard deviation is 1.0.
00788         */
00789     NormalRandomFunctor(Engine const & generator = Engine::global())
00790     : mean_(0.0),
00791       stddev_(1.0),
00792       generator_(generator)
00793     {}
00794 
00795         /** Create functor for normal random numbers with goven mean and standard deviation
00796             using the given engine.
00797         */
00798     NormalRandomFunctor(double mean, double stddev, 
00799                         Engine & generator = Engine::global())
00800     : mean_(mean),
00801       stddev_(stddev),
00802       generator_(generator)
00803     {
00804         vigra_precondition(stddev > 0.0,
00805           "NormalRandomFunctor(): standard deviation must be positive."); 
00806     }
00807     
00808         /** Return a random number as specified in the constructor.
00809         */
00810     double operator()() const
00811     {
00812         return generator_.normal() * stddev_ + mean_;
00813     }
00814 
00815 };
00816 
00817 template <class Engine>
00818 class FunctorTraits<NormalRandomFunctor<Engine> >
00819 {
00820   public:
00821     typedef UniformRandomFunctor<Engine>  type;
00822     
00823     typedef VigraTrueType  isInitializer;
00824     
00825     typedef VigraFalseType isUnaryFunctor;
00826     typedef VigraFalseType isBinaryFunctor;
00827     typedef VigraFalseType isTernaryFunctor;
00828     
00829     typedef VigraFalseType isUnaryAnalyser;
00830     typedef VigraFalseType isBinaryAnalyser;
00831     typedef VigraFalseType isTernaryAnalyser;
00832 };
00833 
00834 //@}
00835 
00836 } // namespace vigra 
00837 
00838 #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.7.1 (3 Dec 2010)