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

vigra/random.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*                  Copyright 2008 by Ullrich Koethe                    */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036
    00037 #ifndef VIGRA_RANDOM_HXX
    00038 #define VIGRA_RANDOM_HXX
    00039
    00040 #include "mathutil.hxx"
    00041 #include "functortraits.hxx"
    00042
    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)