Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkMersenneTwisterRandomVariateGenerator.h

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

Generated at Mon Mar 12 01:37:41 2007 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000