itkMersenneTwisterRandomVariateGenerator.h
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkMersenneTwisterRandomVariateGenerator_h
00018 #define __itkMersenneTwisterRandomVariateGenerator_h
00019
00020 #if defined(_MSC_VER)
00021 #pragma warning ( disable : 4146 )
00022 #endif
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
00112 typedef ITK_UINT32 IntegerType;
00113
00115 itkTypeMacro(MersenneTwisterRandomVariateGenerator,
00116 RandomVariateGeneratorBase );
00117
00119 static Pointer New();
00120 static Pointer GetInstance();
00122
00124 itkStaticConstMacro(StateVectorLength,IntegerType,624);
00125
00127
00128
00129
00130
00131
00132
00134 void Initialize( const IntegerType oneSeed );
00135
00137
00138
00139
00140
00141 void Initialize();
00142
00143
00144
00146 double GetVariateWithClosedRange();
00147
00149 double GetVariateWithClosedRange( const double& n );
00150
00152 double GetVariateWithOpenUpperRange();
00153
00155 double GetVariateWithOpenUpperRange( const double& n );
00156
00158 double GetVariateWithOpenRange();
00159
00161 double GetVariateWithOpenRange( const double& n );
00162
00164 IntegerType GetIntegerVariate();
00165
00167 IntegerType GetIntegerVariate( const IntegerType& n );
00168
00171 double Get53BitVariate();
00172
00173
00174
00175 double GetNormalVariate( const double& mean = 0.0,
00176 const double& variance = 1.0 );
00177
00178
00179
00180 double GetUniformVariate( const double& a, const double& b );
00181
00187 virtual double GetVariate();
00188
00190 double operator()();
00191
00192
00193
00194 inline void SetSeed( const IntegerType oneSeed );
00195 inline void SetSeed( IntegerType * bigSeed, const IntegerType seedLength = StateVectorLength );
00196 inline void SetSeed();
00197
00198
00199
00200
00201
00202
00203
00204 protected:
00205 inline MersenneTwisterRandomVariateGenerator();
00206 virtual ~MersenneTwisterRandomVariateGenerator() {};
00207 virtual void PrintSelf(std::ostream& os, Indent indent) const {
00208 Superclass::PrintSelf(os, indent);
00209
00210
00211 os << indent << "State vector: " << state << std::endl;
00212 os << indent;
00213 register const IntegerType *s = state;
00214 register int i = StateVectorLength;
00215 for(; i--; os << *s++ << "\t" ) {}
00216 os << std::endl;
00217
00218
00219 os << indent << "Next value to be gotten from state: " << pNext << std::endl;
00220
00221
00222 os << indent << "Values left before next reload: " << left << std::endl;
00223 }
00224
00225
00226
00227 itkStaticConstMacro ( M, unsigned int, 397 );
00228
00229 IntegerType state[StateVectorLength];
00230 IntegerType *pNext;
00231 int left;
00232
00233
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 IntegerType twist( const IntegerType& m, const IntegerType& s0, const IntegerType& s1 ) const
00244 {
00245 return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL);
00246 }
00247 static IntegerType hash( vcl_time_t t, vcl_clock_t c );
00248 static Pointer m_Instance;
00249 };
00250
00251
00252
00253
00254 inline MersenneTwisterRandomVariateGenerator::IntegerType
00255 MersenneTwisterRandomVariateGenerator::hash( vcl_time_t t, vcl_clock_t c )
00256 {
00257
00258
00259
00260
00261 static IntegerType differ = 0;
00262
00263 IntegerType h1 = 0;
00264 unsigned char *p = (unsigned char *) &t;
00265 for( size_t i = 0; i < sizeof(t); ++i )
00266 {
00267 h1 *= UCHAR_MAX + 2U;
00268 h1 += p[i];
00269 }
00270 IntegerType h2 = 0;
00271 p = (unsigned char *) &c;
00272 for( size_t j = 0; j < sizeof(c); ++j )
00273 {
00274 h2 *= UCHAR_MAX + 2U;
00275 h2 += p[j];
00276 }
00277 return ( h1 + differ++ ) ^ h2;
00278 }
00279
00280
00281 inline void
00282 MersenneTwisterRandomVariateGenerator::Initialize( const IntegerType seed )
00283 {
00284
00285
00286
00287
00288 register IntegerType *s = state;
00289 register IntegerType *r = state;
00290 register IntegerType i = 1;
00291 *s++ = seed & 0xffffffffUL;
00292 for( i = 1; i < MersenneTwisterRandomVariateGenerator::StateVectorLength; ++i )
00293 {
00294 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
00295 r++;
00296 }
00297 }
00298
00299 inline void
00300 MersenneTwisterRandomVariateGenerator::reload()
00301 {
00302
00303
00304
00305
00306
00307 register int index = static_cast< int >(
00308 M-MersenneTwisterRandomVariateGenerator::StateVectorLength);
00309
00310 register IntegerType *p = state;
00311 register int i;
00312 for( i = MersenneTwisterRandomVariateGenerator::StateVectorLength - M; i--; ++p )
00313 {
00314 *p = twist( p[M], p[0], p[1] );
00315 }
00316 for( i = M; --i; ++p )
00317 {
00318 *p = twist( p[index], p[0], p[1] );
00319 }
00320 *p = twist( p[index], p[0], state[0] );
00321
00322 left = MersenneTwisterRandomVariateGenerator::StateVectorLength, pNext = state;
00323 }
00324
00325 #define SVL 624
00326 inline void
00327 MersenneTwisterRandomVariateGenerator::SetSeed(
00328 IntegerType * const bigSeed, const IntegerType seedLength )
00329 {
00330
00331
00332
00333
00334
00335
00336
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;
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
00387 Initialize(oneSeed);
00388 reload();
00389 }
00390
00391
00392 inline void
00393 MersenneTwisterRandomVariateGenerator::SetSeed()
00394 {
00395
00396 SetSeed( hash( vcl_time(0), vcl_clock() ) );
00397 }
00398
00400 inline MersenneTwisterRandomVariateGenerator::IntegerType
00401 MersenneTwisterRandomVariateGenerator::GetIntegerVariate()
00402 {
00403 if( left == 0 ) reload();
00404 --left;
00406
00407 register IntegerType s1;
00408 s1 = *pNext++;
00409 s1 ^= (s1 >> 11);
00410 s1 ^= (s1 << 7) & 0x9d2c5680UL;
00411 s1 ^= (s1 << 15) & 0xefc60000UL;
00412 return ( s1 ^ (s1 >> 18) );
00413 }
00414
00415
00416 inline double
00417 MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange()
00418 {
00419 return double(GetIntegerVariate()) * (1.0/4294967295.0);
00420 }
00421
00423 inline double
00424 MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange(
00425 const double& n )
00426 {
00427 return GetVariateWithClosedRange() * n;
00428 }
00429
00431 inline double
00432 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange()
00433 {
00434 return double(GetIntegerVariate()) * (1.0/4294967296.0);
00435 }
00436
00438 inline double
00439 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange(
00440 const double& n )
00441 {
00442 return GetVariateWithOpenUpperRange() * n;
00443 }
00444
00446 inline double
00447 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange()
00448 {
00449 return ( double(GetIntegerVariate()) + 0.5 ) * (1.0/4294967296.0);
00450 }
00451
00452
00454 inline double
00455 MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange(
00456 const double& n )
00457 {
00458 return GetVariateWithOpenRange() * n;
00459 }
00460
00461
00462 inline MersenneTwisterRandomVariateGenerator::IntegerType
00463 MersenneTwisterRandomVariateGenerator::GetIntegerVariate(
00464 const IntegerType& n )
00465 {
00466
00467
00468 IntegerType used = n;
00469 used |= used >> 1;
00470 used |= used >> 2;
00471 used |= used >> 4;
00472 used |= used >> 8;
00473 used |= used >> 16;
00474
00475
00476 IntegerType i;
00477 do
00478 {
00479 i = GetIntegerVariate() & used;
00480 }
00481 while( i > n );
00482
00483 return i;
00484 }
00485
00486
00489 inline double
00490 MersenneTwisterRandomVariateGenerator::Get53BitVariate()
00491 {
00492 IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
00493 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
00494 }
00496
00497
00498
00499
00500 inline double
00501 MersenneTwisterRandomVariateGenerator::GetNormalVariate(
00502 const double& mean, const double& variance )
00503 {
00504
00505
00506 double r = vcl_sqrt( -2.0 * vcl_log( 1.0-GetVariateWithOpenRange()) * variance);
00507 double phi = 2.0 * 3.14159265358979323846264338328
00508 * GetVariateWithOpenUpperRange();
00509 return mean + r * vcl_cos(phi);
00510 }
00511
00512
00513
00514
00515 inline double
00516 MersenneTwisterRandomVariateGenerator::GetUniformVariate(
00517 const double& a, const double& b )
00518 {
00519 double u = GetVariateWithOpenUpperRange();
00520 return ((1.0 -u) * a + u * b);
00521 }
00522
00523
00524 inline double
00525 MersenneTwisterRandomVariateGenerator::GetVariate()
00526 {
00527 return GetVariateWithClosedRange();
00528 }
00529
00530
00531 inline double
00532 MersenneTwisterRandomVariateGenerator::operator()()
00533 {
00534 return GetVariate();
00535 }
00536
00537
00538 inline
00539 MersenneTwisterRandomVariateGenerator::MersenneTwisterRandomVariateGenerator()
00540 {
00541 SetSeed ( 121212 );
00542 }
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 }
00587 }
00588
00589 #endif
00590