vigra - vigra/random.hxx Source File

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

vigra/random.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)