ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkMersenneTwisterRandomVariateGenerator.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkMersenneTwisterRandomVariateGenerator_h
00019 #define __itkMersenneTwisterRandomVariateGenerator_h
00020 
00021 #include "itkMacro.h"
00022 #include "itkObjectFactory.h"
00023 #include "itkRandomVariateGeneratorBase.h"
00024 #include "itkIntTypes.h"
00025 #include "vcl_ctime.h"
00026 #include "vnl/vnl_math.h"
00027 #include <limits.h>
00028 
00029 namespace itk
00030 {
00031 namespace Statistics
00032 {
00105 class ITKCommon_EXPORT MersenneTwisterRandomVariateGenerator:
00106   public RandomVariateGeneratorBase
00107 {
00108 public:
00109 
00111   typedef MersenneTwisterRandomVariateGenerator Self;
00112   typedef RandomVariateGeneratorBase            Superclass;
00113   typedef SmartPointer< Self >                  Pointer;
00114   typedef SmartPointer< const Self >            ConstPointer;
00115 
00116   typedef uint32_t IntegerType;
00117 
00119   itkTypeMacro(MersenneTwisterRandomVariateGenerator,
00120                RandomVariateGeneratorBase);
00121 
00125   static Pointer New();
00126 
00132   static Pointer GetInstance();
00133 
00135   itkStaticConstMacro(StateVectorLength, IntegerType, 624);
00136 
00138   void Initialize(const IntegerType oneSeed);
00139 
00140   /* Initialize with vcl_clock time */
00141   void Initialize();
00142 
00144   double GetVariateWithClosedRange();
00145 
00147   double GetVariateWithClosedRange(const double & n);
00148 
00150   double GetVariateWithOpenUpperRange();
00151 
00153   double GetVariateWithOpenUpperRange(const double & n);
00154 
00156   double GetVariateWithOpenRange();
00157 
00159   double GetVariateWithOpenRange(const double & n);
00160 
00162   IntegerType GetIntegerVariate();
00163 
00165   IntegerType GetIntegerVariate(const IntegerType & n);
00166 
00169   double Get53BitVariate();
00170 
00171   /* Access to a normal random number distribution
00172    * TODO: Compare with vnl_sample_normal */
00173   double GetNormalVariate(const double & mean = 0.0,
00174                           const double & variance = 1.0);
00175 
00176   /* Access to a uniform random number distribution in the range [a, b)
00177    * TODO: Compare with vnl_sample_uniform */
00178   double GetUniformVariate(const double & a, const double & b);
00179 
00185   virtual double GetVariate();
00186 
00188   double operator()();
00189 
00190   // Re-seeding functions with same behavior as initializers
00191   inline void SetSeed(const IntegerType oneSeed);
00192   // inline void SetSeed(IntegerType *bigSeed, const IntegerType seedLength = StateVectorLength);
00193   inline void SetSeed();
00194   // Return the current seed
00195   IntegerType GetSeed() { return this->m_Seed; };
00196 
00197   /*
00198   // Saving and loading generator state
00199   void save( IntegerType* saveArray ) const;  // to array of size SAVE
00200   void load( IntegerType *const loadArray );  // from such array
00201   */
00202 protected:
00203   inline MersenneTwisterRandomVariateGenerator();
00204   virtual ~MersenneTwisterRandomVariateGenerator() {}
00205   virtual void PrintSelf(std::ostream & os, Indent indent) const
00206   {
00207     Superclass::PrintSelf(os, indent);
00208 
00209     // Print state vector contents
00210     os << indent << "State vector: " << state << std::endl;
00211     os << indent;
00212     register const IntegerType *s = state;
00213     register int                i = StateVectorLength;
00214     for (; i--; os << *s++ << "\t" ) {}
00215     os << std::endl;
00216 
00217     //Print next value to be gotten from state
00218     os << indent << "Next value to be gotten from state: " << pNext << std::endl;
00219 
00220     //Number of values left before reload
00221     os << indent << "Values left before next reload: " << left << std::endl;
00222   }
00223 
00224   // period parameter
00225   itkStaticConstMacro (M, unsigned int, 397);
00226 
00227   IntegerType  state[StateVectorLength]; // internal state
00228   IntegerType *pNext;                    // next value to get from state
00229   int          left;                     // number of values left before reload
00230                                          // needed
00231   IntegerType  m_Seed;                   // Seed value
00232 
00233   /* Reload array with N new values */
00234   void reload();
00235 
00236   IntegerType hiBit(const IntegerType & u) const { return u & 0x80000000UL; }
00237   IntegerType loBit(const IntegerType & u) const { return u & 0x00000001UL; }
00238   IntegerType loBits(const IntegerType & u) const { return u & 0x7fffffffUL; }
00239   IntegerType mixBits(const IntegerType & u, const IntegerType & v) const
00240   {
00241     return hiBit(u) | loBits(v);
00242   }
00243 
00244   IntegerType twist(const IntegerType & m, const IntegerType & s0, const IntegerType & s1) const
00245   {
00246     return m ^ ( mixBits(s0, s1) >> 1 ) ^ ( -loBit(s1) & 0x9908b0dfUL );
00247   }
00248 
00249   static IntegerType hash(vcl_time_t t, vcl_clock_t c);
00250 
00251   static Pointer m_Instance;
00252 };  // end of class
00253 
00254 // Declare inlined functions.... (must be declared in the header)
00255 
00256 inline MersenneTwisterRandomVariateGenerator::IntegerType
00257 MersenneTwisterRandomVariateGenerator::hash(vcl_time_t t, vcl_clock_t c)
00258 {
00259   // Get a IntegerType from t and c
00260   // Better than IntegerType(x) in case x is floating point in [0,1]
00261   // Based on code by Lawrence Kirby: fred at genesis dot demon dot co dot uk
00262 
00263   static IntegerType differ = 0;  // guarantee time-based seeds will change
00264 
00265   IntegerType    h1 = 0;
00266   unsigned char *p = (unsigned char *)&t;
00267 
00268   const unsigned int sizeOfT = static_cast< unsigned int >( sizeof(t) );
00269   for ( unsigned int i = 0; i < sizeOfT; ++i )
00270     {
00271     h1 *= UCHAR_MAX + 2U;
00272     h1 += p[i];
00273     }
00274   IntegerType h2 = 0;
00275   p = (unsigned char *)&c;
00276 
00277   const unsigned int sizeOfC = static_cast< unsigned int >( sizeof(c) );
00278   for ( unsigned int j = 0; j < sizeOfC; ++j )
00279     {
00280     h2 *= UCHAR_MAX + 2U;
00281     h2 += p[j];
00282     }
00283   return ( h1 + differ++ ) ^ h2;
00284 }
00285 
00286 inline void
00287 MersenneTwisterRandomVariateGenerator::Initialize(const IntegerType seed)
00288 {
00289   this->m_Seed = seed;
00290   // Initialize generator state with seed
00291   // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
00292   // In previous versions, most significant bits (MSBs) of the seed affect
00293   // only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
00294   register IntegerType *s = state;
00295   register IntegerType *r = state;
00296   register IntegerType  i = 1;
00297 
00298   *s++ = seed & 0xffffffffUL;
00299   for ( i = 1; i < MersenneTwisterRandomVariateGenerator::StateVectorLength; ++i )
00300     {
00301     *s++ = ( 1812433253UL * ( *r ^ ( *r >> 30 ) ) + i ) & 0xffffffffUL;
00302     r++;
00303     }
00304 }
00305 
00306 inline void
00307 MersenneTwisterRandomVariateGenerator::reload()
00308 {
00309   // Generate N new values in state
00310   // Made clearer and faster by Matthew Bellew
00311   // matthew dot bellew at home dot com
00312 
00313   // get rid of VS warning
00314   register int index = static_cast< int >(
00315     M - MersenneTwisterRandomVariateGenerator::StateVectorLength );
00316 
00317   register IntegerType *p = state;
00318   register int          i;
00319 
00320   for ( i = MersenneTwisterRandomVariateGenerator::StateVectorLength - M; i--; ++p )
00321     {
00322     *p = twist(p[M], p[0], p[1]);
00323     }
00324   for ( i = M; --i; ++p )
00325     {
00326     *p = twist(p[index], p[0], p[1]);
00327     }
00328   *p = twist(p[index], p[0], state[0]);
00329 
00330   left = MersenneTwisterRandomVariateGenerator::StateVectorLength, pNext = state;
00331 }
00332 
00333   /*
00334 #define SVL 624
00335 inline void
00336 MersenneTwisterRandomVariateGenerator::SetSeed(
00337   IntegerType *const bigSeed, const IntegerType seedLength)
00338 {
00339   // Seed the generator with an array of IntegerType's
00340   // There are 2^19937-1 possible initial states.  This function allows
00341   // all of those to be accessed by providing at least 19937 bits (with a
00342   // default seed length of StateVectorLength = 624 IntegerType's).
00343   // Any bits above the lower 32
00344   // in each element are discarded.
00345   // Just call seed() if you want to get array from /dev/urandom
00346   Initialize(19650218UL);
00347   register IntegerType i = 1;
00348   register IntegerType j = 0;
00349   register int         k;
00350   if ( StateVectorLength > seedLength )
00351     {
00352     k = StateVectorLength;
00353     }
00354   else
00355     {
00356     k = seedLength;
00357     }
00358   for (; k; --k )
00359     {
00360     state[i] =
00361       state[i] ^ ( ( state[i - 1] ^ ( state[i - 1] >> 30 ) ) * 1664525UL );
00362     state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
00363     state[i] &= 0xffffffffUL;
00364     ++i;  ++j;
00365     if ( i >= StateVectorLength ) { state[0] = state[StateVectorLength - 1];  i = 1; }
00366     if ( j >= seedLength ) { j = 0; }
00367     }
00368   for ( k = StateVectorLength - 1; k; --k )
00369     {
00370     state[i] =
00371       state[i] ^ ( ( state[i - 1] ^ ( state[i - 1] >> 30 ) ) * 1566083941UL );
00372     state[i] -= i;
00373     state[i] &= 0xffffffffUL;
00374     ++i;
00375     if ( i >= SVL )
00376       {
00377       state[0] = state[StateVectorLength - 1];  i = 1;
00378       }
00379     }
00380   state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
00381   reload();
00382 }
00383   */
00384 
00385 inline void
00386 MersenneTwisterRandomVariateGenerator::Initialize()
00387 {
00388   SetSeed();
00389 }
00390 
00391 inline void
00392 MersenneTwisterRandomVariateGenerator::SetSeed(const IntegerType oneSeed)
00393 {
00394   // Seed the generator with a simple IntegerType
00395   Initialize(oneSeed);
00396   reload();
00397 }
00398 
00399 inline void
00400 MersenneTwisterRandomVariateGenerator::SetSeed()
00401 {
00402   // use time() and clock() to generate a unlikely-to-repeat seed.
00403   SetSeed( hash( vcl_time(0), vcl_clock() ) );
00404 }
00405 
00407 inline MersenneTwisterRandomVariateGenerator::IntegerType
00408 MersenneTwisterRandomVariateGenerator::GetIntegerVariate()
00409 {
00410   if ( left == 0 ) { reload(); }
00411   --left;
00413 
00414   register IntegerType s1;
00415   s1 = *pNext++;
00416   s1 ^= ( s1 >> 11 );
00417   s1 ^= ( s1 <<  7 ) & 0x9d2c5680UL;
00418   s1 ^= ( s1 << 15 ) & 0xefc60000UL;
00419   return ( s1 ^ ( s1 >> 18 ) );
00420 }
00421 
00422 inline double
00423 MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange()
00424 {
00425   return double( GetIntegerVariate() ) * ( 1.0 / 4294967295.0 );
00426 }
00427 
00429 inline double
00430 MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange(
00431   const double & n)
00432 {
00433   return GetVariateWithClosedRange() * n;
00434 }
00435 
00437 inline double
00438 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange()
00439 {
00440   return double( GetIntegerVariate() ) * ( 1.0 / 4294967296.0 );
00441 }
00442 
00444 inline double
00445 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange(
00446   const double & n)
00447 {
00448   return GetVariateWithOpenUpperRange() * n;
00449 }
00450 
00452 inline double
00453 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange()
00454 {
00455   return ( double( GetIntegerVariate() ) + 0.5 ) * ( 1.0 / 4294967296.0 );
00456 }
00457 
00459 inline double
00460 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange(
00461   const double & n)
00462 {
00463   return GetVariateWithOpenRange() * n;
00464 }
00465 
00466 inline MersenneTwisterRandomVariateGenerator::IntegerType
00467 MersenneTwisterRandomVariateGenerator::GetIntegerVariate(
00468   const IntegerType & n)
00469 {
00470   // Find which bits are used in n
00471   // Optimized by Magnus Jonsson magnus at smartelectronix dot com
00472   IntegerType used = n;
00473 
00474   used |= used >> 1;
00475   used |= used >> 2;
00476   used |= used >> 4;
00477   used |= used >> 8;
00478   used |= used >> 16;
00479 
00480   // Draw numbers until one is found in [0,n]
00481   IntegerType i;
00482   do
00483     {
00484     i = GetIntegerVariate() & used;  // toss unused bits to shorten search
00485     }
00486   while ( i > n );
00487 
00488   return i;
00489 }
00490 
00493 inline double
00494 MersenneTwisterRandomVariateGenerator::Get53BitVariate()
00495 {
00496   IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
00497 
00498   return ( a * 67108864.0 + b ) * ( 1.0 / 9007199254740992.0 );  // by Isaku
00499                                                                  // Wada
00500 }
00501 
00502 /* Access to a normal random number distribution */
00503 // TODO: Compare with vnl_sample_normal
00504 inline double
00505 MersenneTwisterRandomVariateGenerator::GetNormalVariate(
00506   const double & mean, const double & variance)
00507 {
00508   // Return a real number from a normal (Gaussian) distribution with given
00509   // mean and variance by Box-Muller method
00510   double r = vcl_sqrt(-2.0 * vcl_log( 1.0 - GetVariateWithOpenRange() ) * variance);
00511   double phi = 2.0 * vnl_math::pi
00512                * GetVariateWithOpenUpperRange();
00513 
00514   return mean + r *vcl_cos(phi);
00515 }
00516 
00517 /* Access to a uniform random number distribution */
00518 // TODO: Compare with vnl_sample_uniform
00519 inline double
00520 MersenneTwisterRandomVariateGenerator::GetUniformVariate(
00521   const double & a, const double & b)
00522 {
00523   double u = GetVariateWithOpenUpperRange();
00524 
00525   return ( ( 1.0 - u ) * a + u * b );
00526 }
00527 
00528 inline double
00529 MersenneTwisterRandomVariateGenerator::GetVariate()
00530 {
00531   return GetVariateWithClosedRange();
00532 }
00533 
00534 inline double
00535 MersenneTwisterRandomVariateGenerator::operator()()
00536 {
00537   return GetVariate();
00538 }
00539 
00540 inline
00541 MersenneTwisterRandomVariateGenerator::MersenneTwisterRandomVariateGenerator()
00542 {
00543   SetSeed (121212);
00544 }
00545 
00546 /* Change log from MTRand.h */
00547 // Change log:
00548 //
00549 // v0.1 - First release on 15 May 2000
00550 //      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
00551 //      - Translated from C to C++
00552 //      - Made completely ANSI compliant
00553 //      - Designed convenient interface for initialization, seeding, and
00554 //        obtaining numbers in default or user-defined ranges
00555 //      - Added automatic seeding from /dev/urandom or time() and clock()
00556 //      - Provided functions for saving and loading generator state
00557 //
00558 // v0.2 - Fixed bug which reloaded generator one step too late
00559 //
00560 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
00561 //
00562 // v0.4 - Removed trailing newline in saved generator format to be consistent
00563 //        with output format of built-in types
00564 //
00565 // v0.5 - Improved portability by replacing static const int's with enum's and
00566 //        clarifying return values in seed(); suggested by Eric Heimburg
00567 //      - Removed MAXINT constant; use 0xffffffffUL instead
00568 //
00569 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
00570 //      - Changed integer [0,n] generator to give better uniformity
00571 //
00572 // v0.7 - Fixed operator precedence ambiguity in reload()
00573 //      - Added access for real numbers in (0,1) and (0,n)
00574 //
00575 // v0.8 - Included time.h header to properly support time_t and clock_t
00576 //
00577 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
00578 //      - Allowed for seeding with arrays of any length
00579 //      - Added access for real numbers in [0,1) with 53-bit resolution
00580 //      - Added access for real numbers from normal (Gaussian) distributions
00581 //      - Increased overall speed by optimizing twist()
00582 //      - Doubled speed of integer [0,n] generation
00583 //      - Fixed out-of-range number generation on 64-bit machines
00584 //      - Improved portability by substituting literal constants for long enum's
00585 //      - Changed license from GNU LGPL to BSD
00586 } // end of namespace Statistics
00587 } // end of namespace itk
00588 
00589 #endif
00590