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
00103 class ITKCommon_EXPORT MersenneTwisterRandomVariateGenerator :
00104 public RandomVariateGeneratorBase
00105 {
00106 public:
00107
00109 typedef MersenneTwisterRandomVariateGenerator Self;
00110 typedef RandomVariateGeneratorBase Superclass;
00111 typedef SmartPointer<Self> Pointer;
00112 typedef SmartPointer<const Self> ConstPointer;
00113
00114 typedef uint32_t IntegerType;
00115
00117 itkTypeMacro(MersenneTwisterRandomVariateGenerator,
00118 RandomVariateGeneratorBase );
00119
00121 static Pointer New();
00122 static Pointer GetInstance();
00124
00126 itkStaticConstMacro(StateVectorLength,IntegerType,624);
00127
00129
00130
00131
00132
00133
00134
00136 void Initialize( const IntegerType oneSeed );
00137
00139
00140
00141
00142
00143 void Initialize();
00144
00145
00146
00148 double GetVariateWithClosedRange();
00149
00151 double GetVariateWithClosedRange( const double& n );
00152
00154 double GetVariateWithOpenUpperRange();
00155
00157 double GetVariateWithOpenUpperRange( const double& n );
00158
00160 double GetVariateWithOpenRange();
00161
00163 double GetVariateWithOpenRange( const double& n );
00164
00166 IntegerType GetIntegerVariate();
00167
00169 IntegerType GetIntegerVariate( const IntegerType& n );
00170
00173 double Get53BitVariate();
00174
00175
00176
00177 double GetNormalVariate( const double& mean = 0.0,
00178 const double& variance = 1.0 );
00179
00180
00181
00182 double GetUniformVariate( const double& a, const double& b );
00183
00189 virtual double GetVariate();
00190
00192 double operator()();
00193
00194
00195
00196 inline void SetSeed( const IntegerType oneSeed );
00197 inline void SetSeed( IntegerType * bigSeed, const IntegerType seedLength = StateVectorLength );
00198 inline void SetSeed();
00199
00200
00201
00202
00203
00204
00205
00206 protected:
00207 inline MersenneTwisterRandomVariateGenerator();
00208 virtual ~MersenneTwisterRandomVariateGenerator() {};
00209 virtual void PrintSelf(std::ostream& os, Indent indent) const {
00210 Superclass::PrintSelf(os, indent);
00211
00212
00213 os << indent << "State vector: " << state << std::endl;
00214 os << indent;
00215 register const IntegerType *s = state;
00216 register int i = StateVectorLength;
00217 for(; i--; os << *s++ << "\t" ) {}
00218 os << std::endl;
00219
00220
00221 os << indent << "Next value to be gotten from state: " << pNext << std::endl;
00222
00223
00224 os << indent << "Values left before next reload: " << left << std::endl;
00225 }
00226
00227
00228
00229 itkStaticConstMacro ( M, unsigned int, 397 );
00230
00231 IntegerType state[StateVectorLength];
00232 IntegerType *pNext;
00233 int left;
00234
00235
00236 void reload();
00237
00238 IntegerType hiBit( const IntegerType& u ) const { return u & 0x80000000UL; }
00239 IntegerType loBit( const IntegerType& u ) const { return u & 0x00000001UL; }
00240 IntegerType loBits( const IntegerType& u ) const { return u & 0x7fffffffUL; }
00241 IntegerType mixBits( const IntegerType& u, const IntegerType& v ) const
00242 {
00243 return hiBit(u) | loBits(v);
00244 }
00245 IntegerType twist( const IntegerType& m, const IntegerType& s0, const IntegerType& s1 ) const
00246 {
00247 return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL);
00248 }
00249 static IntegerType hash( vcl_time_t t, vcl_clock_t c );
00250 static Pointer m_Instance;
00251 };
00252
00253
00254
00255
00256 inline MersenneTwisterRandomVariateGenerator::IntegerType
00257 MersenneTwisterRandomVariateGenerator::hash( vcl_time_t t, vcl_clock_t c )
00258 {
00259
00260
00261
00262
00263 static IntegerType differ = 0;
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
00287
00288
00289
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
00305
00306
00307
00308
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 #define SVL 624
00328 inline void
00329 MersenneTwisterRandomVariateGenerator::SetSeed(
00330 IntegerType * const bigSeed, const IntegerType seedLength )
00331 {
00332
00333
00334
00335
00336
00337
00338
00339 Initialize(19650218UL);
00340 register IntegerType i = 1;
00341 register IntegerType j = 0;
00342 register int k;
00343 if ( StateVectorLength > seedLength )
00344 {
00345 k = StateVectorLength;
00346 }
00347 else
00348 {
00349 k = seedLength;
00350 }
00351 for(; k; --k )
00352 {
00353 state[i] =
00354 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
00355 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
00356 state[i] &= 0xffffffffUL;
00357 ++i; ++j;
00358 if( i >= StateVectorLength ) { state[0] = state[StateVectorLength-1]; i = 1; }
00359 if( j >= seedLength ) j = 0;
00360 }
00361 for( k = StateVectorLength - 1; k; --k )
00362 {
00363 state[i] =
00364 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
00365 state[i] -= i;
00366 state[i] &= 0xffffffffUL;
00367 ++i;
00368 if( i >= SVL )
00369 {
00370 state[0] = state[StateVectorLength-1]; i = 1;
00371 }
00372 }
00373 state[0] = 0x80000000UL;
00374 reload();
00375 }
00376
00377
00378 inline void
00379 MersenneTwisterRandomVariateGenerator::Initialize()
00380 {
00381 SetSeed();
00382 }
00383
00384
00385 inline void
00386 MersenneTwisterRandomVariateGenerator::SetSeed( const IntegerType oneSeed )
00387 {
00388
00389 Initialize(oneSeed);
00390 reload();
00391 }
00392
00393
00394 inline void
00395 MersenneTwisterRandomVariateGenerator::SetSeed()
00396 {
00397
00398 SetSeed( hash( vcl_time(0), vcl_clock() ) );
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
00469
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
00478 IntegerType i;
00479 do
00480 {
00481 i = GetIntegerVariate() & used;
00482 }
00483 while( i > n );
00484
00485 return i;
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);
00496 }
00498
00499
00500
00501
00502 inline double
00503 MersenneTwisterRandomVariateGenerator::GetNormalVariate(
00504 const double& mean, const double& variance )
00505 {
00506
00507
00508 double r = vcl_sqrt( -2.0 * vcl_log( 1.0-GetVariateWithOpenRange()) * variance);
00509 double phi = 2.0 * vnl_math::pi
00510 * GetVariateWithOpenUpperRange();
00511 return mean + r * vcl_cos(phi);
00512 }
00513
00514
00515
00516
00517 inline double
00518 MersenneTwisterRandomVariateGenerator::GetUniformVariate(
00519 const double& a, const double& b )
00520 {
00521 double u = GetVariateWithOpenUpperRange();
00522 return ((1.0 -u) * a + u * b);
00523 }
00524
00525
00526 inline double
00527 MersenneTwisterRandomVariateGenerator::GetVariate()
00528 {
00529 return GetVariateWithClosedRange();
00530 }
00531
00532
00533 inline double
00534 MersenneTwisterRandomVariateGenerator::operator()()
00535 {
00536 return GetVariate();
00537 }
00538
00539
00540 inline
00541 MersenneTwisterRandomVariateGenerator::MersenneTwisterRandomVariateGenerator()
00542 {
00543 SetSeed ( 121212 );
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 }
00590
00591 #endif
00592