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: 2008-04-07 12:18:26 $
00007   Version:   $Revision: 1.6 $
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 #if defined(_MSC_VER)
00018 #pragma warning ( disable : 4146 )
00019 #endif
00020 
00021 #ifndef __itkMersenneTwisterRandomVariateGenerator_h
00022 #define __itkMersenneTwisterRandomVariateGenerator_h
00023 
00024 #include "itkMacro.h"
00025 #include "itkObjectFactory.h"
00026 #include "itkRandomVariateGeneratorBase.h"
00027 #include "itkIntTypes.h"
00028 #include "vcl_ctime.h"
00029 #include "vnl/vnl_math.h"
00030 #include <limits.h>
00031 
00032 
00033 namespace itk {
00034 namespace Statistics {
00035 
00101 class ITKCommon_EXPORT MersenneTwisterRandomVariateGenerator : 
00102     public RandomVariateGeneratorBase
00103 {
00104 public:
00105   
00107   typedef MersenneTwisterRandomVariateGenerator Self ;
00108   typedef RandomVariateGeneratorBase Superclass;
00109   typedef SmartPointer<Self>   Pointer;
00110   typedef SmartPointer<const Self>  ConstPointer;
00111   typedef ITK_UINT32 IntegerType;
00112 
00114   itkTypeMacro(MersenneTwisterRandomVariateGenerator, 
00115                RandomVariateGeneratorBase );
00116 
00118   static Pointer New();
00119   static Pointer GetInstance();
00121 
00123   itkStaticConstMacro(StateVectorLength,IntegerType,624);
00124 
00126   // itkStaticConstMacro(SAVE,IntegerType,625);
00127 
00128   
00129   
00130   // Methods to reseed
00131   
00133   void Initialize( const IntegerType oneSeed );
00134 
00136   //void Initialize( IntegerType *const bigSeed, 
00137   //        IntegerType const seedLength = StateVectorLength );  
00138 
00139   /* Initialize with vcl_clock time */ 
00140   void Initialize();  
00141   
00142   // Methods to get various random variates ...
00143   
00145   double GetVariateWithClosedRange();
00146 
00148   double GetVariateWithClosedRange( const double& n );
00149 
00151   double GetVariateWithOpenUpperRange();
00152 
00154   double GetVariateWithOpenUpperRange( const double& n );
00155 
00157   double GetVariateWithOpenRange();
00158 
00160   double GetVariateWithOpenRange( const double& n );
00161 
00163   IntegerType GetIntegerVariate();
00164 
00166   IntegerType GetIntegerVariate( const IntegerType& n );      
00167 
00170   double Get53BitVariate();
00171 
00172   /* Access to a normal random number distribution 
00173    * TODO: Compare with vnl_sample_normal */
00174   double GetNormalVariate( const double& mean = 0.0, 
00175       const double& variance = 1.0 );
00176   
00177   /* Access to a uniform random number distribution in the range [a, b)
00178    * TODO: Compare with vnl_sample_uniform */
00179   double GetUniformVariate( const double& a, const double& b );
00180   
00186   virtual double GetVariate();
00187 
00189   double operator()(); 
00190 
00191 
00192   // Re-seeding functions with same behavior as initializers
00193   inline void SetSeed( const IntegerType oneSeed );
00194   inline void SetSeed( IntegerType * bigSeed, const IntegerType seedLength = StateVectorLength );
00195   inline void SetSeed();
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 
00203  protected:
00204   inline MersenneTwisterRandomVariateGenerator();
00205   virtual ~MersenneTwisterRandomVariateGenerator() {}; 
00206   virtual void PrintSelf(std::ostream& os, Indent indent) const {
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 
00225   // enum { M = 397 };  // period parameter
00226   itkStaticConstMacro ( M, unsigned int, 397 );
00227   
00228   IntegerType state[StateVectorLength];   // internal state
00229   IntegerType *pNext;     // next value to get from state
00230   int left;          // number of values left before reload needed
00231 
00232   /* Reload array with N new values */
00233   void reload();
00234   
00235   IntegerType hiBit( const IntegerType& u ) const { return u & 0x80000000UL; }
00236   IntegerType loBit( const IntegerType& u ) const { return u & 0x00000001UL; }
00237   IntegerType loBits( const IntegerType& u ) const { return u & 0x7fffffffUL; }
00238   IntegerType mixBits( const IntegerType& u, const IntegerType& v ) const
00239     { 
00240     return hiBit(u) | loBits(v); 
00241     }
00242   IntegerType twist( const IntegerType& m, const IntegerType& s0, const IntegerType& s1 ) const
00243     { 
00244     return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); 
00245     }
00246   static IntegerType hash( vcl_time_t t, vcl_clock_t c );
00247   static Pointer m_Instance;
00248 } ;  // end of class
00249   
00250 
00251 
00252 
00253 // Declare inlined functions.... (must be declared in the header)
00254 // Declare then in order to keep SGI happy
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   for( size_t i = 0; i < sizeof(t); ++i )
00268     {
00269     h1 *= UCHAR_MAX + 2U;
00270     h1 += p[i];
00271     }
00272   IntegerType h2 = 0;
00273   p = (unsigned char *) &c;
00274   for( size_t j = 0; j < sizeof(c); ++j )
00275     {
00276     h2 *= UCHAR_MAX + 2U;
00277     h2 += p[j];
00278     }
00279   return ( h1 + differ++ ) ^ h2;
00280   }
00281 
00282 
00283 inline void 
00284   MersenneTwisterRandomVariateGenerator::Initialize( const IntegerType seed )
00285   {
00286   // Initialize generator state with seed
00287   // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
00288   // In previous versions, most significant bits (MSBs) of the seed affect
00289   // only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
00290   register IntegerType *s = state;
00291   register IntegerType *r = state;
00292   register IntegerType i = 1;
00293   *s++ = seed & 0xffffffffUL;
00294   for( i = 1; i < MersenneTwisterRandomVariateGenerator::StateVectorLength; ++i )
00295     {
00296     *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
00297     r++;
00298     }
00299   }
00300 
00301 inline void 
00302   MersenneTwisterRandomVariateGenerator::reload()
00303   {
00304   // Generate N new values in state
00305   // Made clearer and faster by Matthew Bellew 
00306   // matthew dot bellew at home dot com
00307   
00308   // get rid of VS warning
00309   register int index = static_cast< int >(
00310       M-MersenneTwisterRandomVariateGenerator::StateVectorLength);
00311     
00312   register IntegerType *p = state;
00313   register int i;
00314   for( i = MersenneTwisterRandomVariateGenerator::StateVectorLength - M; i--; ++p )
00315     {
00316     *p = twist( p[M], p[0], p[1] );
00317     }
00318   for( i = M; --i; ++p )
00319     {
00320     *p = twist( p[index], p[0], p[1] );
00321     }
00322   *p = twist( p[index], p[0], state[0] );
00323 
00324   left = MersenneTwisterRandomVariateGenerator::StateVectorLength, pNext = state;
00325   }
00326 
00327 
00328 
00329 #define SVL 624
00330 inline void 
00331   MersenneTwisterRandomVariateGenerator::SetSeed( 
00332       IntegerType * const bigSeed, const IntegerType seedLength )
00333   {
00334   // Seed the generator with an array of IntegerType's
00335   // There are 2^19937-1 possible initial states.  This function allows
00336   // all of those to be accessed by providing at least 19937 bits (with a
00337   // default seed length of StateVectorLength = 624 IntegerType's).  
00338   // Any bits above the lower 32
00339   // in each element are discarded.
00340   // Just call seed() if you want to get array from /dev/urandom
00341   Initialize(19650218UL);
00342   register IntegerType i = 1;
00343   register IntegerType j = 0;
00344   register int k;
00345   if ( StateVectorLength > seedLength )
00346     {
00347     k = StateVectorLength;
00348     }
00349   else
00350     {
00351     k = seedLength;
00352     }
00353   for( ; k; --k )
00354     {
00355     state[i] =
00356       state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
00357     state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
00358     state[i] &= 0xffffffffUL;
00359     ++i;  ++j;
00360     if( i >= StateVectorLength ) { state[0] = state[StateVectorLength-1];  i = 1; }
00361     if( j >= seedLength ) j = 0;
00362     }
00363   for( k = StateVectorLength - 1; k; --k )
00364     {
00365     state[i] =
00366       state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
00367     state[i] -= i;
00368     state[i] &= 0xffffffffUL;
00369     ++i;
00370     if( i >= SVL ) 
00371       { 
00372       state[0] = state[StateVectorLength-1];  i = 1; 
00373       }
00374     }
00375   state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
00376   reload();
00377 }
00378 
00379 
00380 inline void
00381   MersenneTwisterRandomVariateGenerator::Initialize()
00382   { 
00383   SetSeed(); 
00384   }
00385 
00386 
00387 inline void 
00388   MersenneTwisterRandomVariateGenerator::SetSeed( const IntegerType oneSeed )
00389   {
00390   // Seed the generator with a simple IntegerType
00391   Initialize(oneSeed);
00392   reload();
00393   }
00394 
00395 
00396 inline void 
00397   MersenneTwisterRandomVariateGenerator::SetSeed()
00398   {
00399   // use time() and clock() to generate a unlikely-to-repeat seed.
00400   SetSeed( hash( vcl_time(0), vcl_clock() ) );
00401   }
00402 
00403 
00404 
00406 inline MersenneTwisterRandomVariateGenerator::IntegerType 
00407   MersenneTwisterRandomVariateGenerator::GetIntegerVariate()
00408   {
00409   if( left == 0 ) reload();
00410   --left;
00412 
00413   register IntegerType s1;
00414   s1 = *pNext++;
00415   s1 ^= (s1 >> 11);
00416   s1 ^= (s1 <<  7) & 0x9d2c5680UL;
00417   s1 ^= (s1 << 15) & 0xefc60000UL;
00418   return ( s1 ^ (s1 >> 18) );
00419   }
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 
00458 
00460 inline double 
00461   MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange( 
00462                                                   const double& n )
00463   { 
00464   return GetVariateWithOpenRange() * n; 
00465   }
00466 
00467 
00468 inline MersenneTwisterRandomVariateGenerator::IntegerType 
00469   MersenneTwisterRandomVariateGenerator::GetIntegerVariate( 
00470                                         const IntegerType& n )
00471   {
00472   // Find which bits are used in n
00473   // Optimized by Magnus Jonsson magnus at smartelectronix dot com
00474   IntegerType used = n;
00475   used |= used >> 1;
00476   used |= used >> 2;
00477   used |= used >> 4;
00478   used |= used >> 8;
00479   used |= used >> 16;
00480 
00481   // Draw numbers until one is found in [0,n]
00482   IntegerType i;
00483   do
00484     {
00485     i = GetIntegerVariate() & used;  // toss unused bits to shorten search
00486     }  while( i > n );
00487   
00488   return i;
00489   }
00490 
00491 
00492 
00495 inline double 
00496   MersenneTwisterRandomVariateGenerator::Get53BitVariate()  
00497   {
00498   IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
00499   return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
00500   }
00502 
00503 
00504 /* Access to a normal random number distribution */
00505 // TODO: Compare with vnl_sample_normal
00506 inline double 
00507   MersenneTwisterRandomVariateGenerator::GetNormalVariate( 
00508       const double& mean, const double& variance )
00509   {
00510   // Return a real number from a normal (Gaussian) distribution with given
00511   // mean and variance by Box-Muller method
00512   double r = vcl_sqrt( -2.0 * vcl_log( 1.0-GetVariateWithOpenRange()) * variance);
00513   double phi = 2.0 * 3.14159265358979323846264338328 
00514                           * GetVariateWithOpenUpperRange();
00515   return mean + r * vcl_cos(phi);
00516   }
00517 
00518 
00519 
00520 /* Access to a uniform random number distribution */
00521 // TODO: Compare with vnl_sample_uniform
00522 inline double 
00523   MersenneTwisterRandomVariateGenerator::GetUniformVariate( 
00524                             const double& a, const double& b )
00525   {
00526   double u = GetVariateWithOpenUpperRange();
00527   return ((1.0 -u) * a + u * b); 
00528   }
00529 
00530 
00531 inline double 
00532   MersenneTwisterRandomVariateGenerator::GetVariate() 
00533   {
00534   return GetVariateWithClosedRange();
00535   }
00536 
00537 
00538 inline double 
00539   MersenneTwisterRandomVariateGenerator::operator()()
00540   { 
00541   return GetVariate(); 
00542   }  
00543  
00544 
00545 inline 
00546   MersenneTwisterRandomVariateGenerator::MersenneTwisterRandomVariateGenerator()
00547   {
00548     SetSeed ( 121212 );
00549   }
00550   
00551 
00552 /* Change log from MTRand.h */
00553 // Change log:
00554 //
00555 // v0.1 - First release on 15 May 2000
00556 //      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
00557 //      - Translated from C to C++
00558 //      - Made completely ANSI compliant
00559 //      - Designed convenient interface for initialization, seeding, and
00560 //        obtaining numbers in default or user-defined ranges
00561 //      - Added automatic seeding from /dev/urandom or time() and clock()
00562 //      - Provided functions for saving and loading generator state
00563 //
00564 // v0.2 - Fixed bug which reloaded generator one step too late
00565 //
00566 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
00567 //
00568 // v0.4 - Removed trailing newline in saved generator format to be consistent
00569 //        with output format of built-in types
00570 //
00571 // v0.5 - Improved portability by replacing static const int's with enum's and
00572 //        clarifying return values in seed(); suggested by Eric Heimburg
00573 //      - Removed MAXINT constant; use 0xffffffffUL instead
00574 //
00575 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
00576 //      - Changed integer [0,n] generator to give better uniformity
00577 //
00578 // v0.7 - Fixed operator precedence ambiguity in reload()
00579 //      - Added access for real numbers in (0,1) and (0,n)
00580 //
00581 // v0.8 - Included time.h header to properly support time_t and clock_t
00582 //
00583 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
00584 //      - Allowed for seeding with arrays of any length
00585 //      - Added access for real numbers in [0,1) with 53-bit resolution
00586 //      - Added access for real numbers from normal (Gaussian) distributions
00587 //      - Increased overall speed by optimizing twist()
00588 //      - Doubled speed of integer [0,n] generation
00589 //      - Fixed out-of-range number generation on 64-bit machines
00590 //      - Improved portability by substituting literal constants for long enum's
00591 //      - Changed license from GNU LGPL to BSD
00592 
00593 } // end of namespace Statistics
00594 } // end of namespace itk
00595 
00596 #endif
00597 
00598 
00599 

Generated at Tue Jul 29 21:14:15 2008 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000