ITK  5.1.0
Insight Toolkit
itkMersenneTwisterRandomVariateGenerator.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkMersenneTwisterRandomVariateGenerator_h
19 #define itkMersenneTwisterRandomVariateGenerator_h
20 
21 #include "itkMacro.h"
22 #include "itkObjectFactory.h"
24 #include "itkIntTypes.h"
25 #include "itkMath.h"
26 #include "itkSingletonMacro.h"
27 
28 #include <atomic>
29 #include <mutex>
30 #include <climits>
31 #include <ctime>
32 
33 namespace itk
34 {
35 namespace Statistics
36 {
125 struct MersenneTwisterGlobals;
126 
128 {
129 public:
135 
137 
140 
149  static Pointer
150  New();
151 
160  static Pointer
161  GetInstance();
162 
164  static constexpr IntegerType StateVectorLength = 624;
165 
167  void
168  Initialize(const IntegerType oneSeed);
169 
170  /* Initialize with clock time */
171  void
172  Initialize();
173 
175  double
176  GetVariateWithClosedRange();
177 
179  double
180  GetVariateWithClosedRange(const double & n);
181 
183  double
184  GetVariateWithOpenUpperRange();
185 
187  double
188  GetVariateWithOpenUpperRange(const double & n);
189 
191  double
192  GetVariateWithOpenRange();
193 
195  double
196  GetVariateWithOpenRange(const double & n);
197 
200  GetIntegerVariate();
201 
204  GetIntegerVariate(const IntegerType & n);
205 
208  double
209  Get53BitVariate();
210 
211  /* Access to a normal random number distribution
212  * TODO: Compare with vnl_sample_normal */
213  double
214  GetNormalVariate(const double & mean = 0.0, const double & variance = 1.0);
215 
216  /* Access to a uniform random number distribution in the range [a, b)
217  * TODO: Compare with vnl_sample_uniform */
218  double
219  GetUniformVariate(const double & a, const double & b);
220 
226  double
227  GetVariate() override;
228 
230  double
231  operator()();
232 
237  inline void
238  SetSeed(const IntegerType oneSeed);
239  inline void
240  SetSeed();
242 
248  GetSeed() const;
249 
255  static IntegerType
256  GetNextSeed();
257 
258  /*
259  // Saving and loading generator state
260  void save( IntegerType* saveArray ) const; // to array of size SAVE
261  void load( IntegerType *const loadArray ); // from such array
262  */
263 
264 protected:
267  void
268  PrintSelf(std::ostream & os, Indent indent) const override;
269 
271  static constexpr unsigned int M = 397;
272 
274  void
275  reload();
276 
278  hiBit(const IntegerType & u) const
279  {
280  return u & 0x80000000;
281  }
282  IntegerType
283  loBit(const IntegerType & u) const
284  {
285  return u & 0x00000001;
286  }
287  IntegerType
288  loBits(const IntegerType & u) const
289  {
290  return u & 0x7fffffff;
291  }
292  IntegerType
293  mixBits(const IntegerType & u, const IntegerType & v) const
294  {
295  return hiBit(u) | loBits(v);
296  }
297 
298  IntegerType
299  twist(const IntegerType & m, const IntegerType & s0, const IntegerType & s1) const
300  {
301  return m ^ (mixBits(s0, s1) >> 1) ^ (-static_cast<int32_t>(loBit(s1)) & 0x9908b0df);
302  }
303 
304  static IntegerType
305  hash(time_t t, clock_t c);
306 
307  // Internal state
308  IntegerType state[StateVectorLength];
309 
310  // Next value to get from state
312 
313  // Number of values left before reload is needed
314  int m_Left;
315 
316  // Seed value
317  std::atomic<IntegerType> m_Seed;
318 
319 private:
321  itkGetGlobalDeclarationMacro(MersenneTwisterGlobals, PimplGlobals);
322 
324  static Pointer
325  CreateInstance();
326 
327  // Local lock to enable concurrent access to singleton
328  std::mutex m_InstanceLock;
329 
330  // Static/Global Variable need to be thread-safely accessed
331 
332  static MersenneTwisterGlobals * m_PimplGlobals;
333 
334 }; // end of class
335 
336 // Declare inlined functions.... (must be declared in the header)
337 
338 inline void
340 {
341  std::lock_guard<std::mutex> mutexHolder(m_InstanceLock);
342  this->m_Seed = seed;
343  // Initialize generator state with seed
344  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
345  // In previous versions, most significant bits (MSBs) of the seed affect
346  // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
347  IntegerType * s = state;
348  IntegerType * r = state;
349  IntegerType i = 1;
350 
351  *s++ = seed & 0xffffffffUL;
353  {
354  *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
355  r++;
356  }
357  reload();
358 }
359 
360 inline void
362 {
363  // Generate N new values in state
364  // Made clearer and faster by Matthew Bellew
365  // matthew dot bellew at home dot com
366 
367  // get rid of VS warning
368  auto index = static_cast<int>(M - MersenneTwisterRandomVariateGenerator::StateVectorLength);
369 
370  IntegerType * p = state;
371  int i;
372 
374  {
375  *p = twist(p[M], p[0], p[1]);
376  }
377  for (i = M; --i; ++p)
378  {
379  *p = twist(p[index], p[0], p[1]);
380  }
381  *p = twist(p[index], p[0], state[0]);
382 
384  m_PNext = state;
385 }
386 
387 inline void
389 {
390  SetSeed();
391 }
392 
393 inline void
395 {
396  // Seed the generator with a simple IntegerType
397  Initialize(oneSeed);
398 }
399 
400 inline void
402 {
403  // use time() and clock() to generate a unlikely-to-repeat seed.
404  SetSeed(hash(time(nullptr), clock()));
405 }
406 
407 
410 {
411  return this->m_Seed;
412 }
413 
417 {
418  if (m_Left == 0)
419  {
420  reload();
421  }
422  --m_Left;
424 
425  IntegerType s1 = *m_PNext++;
426  s1 ^= (s1 >> 11);
427  s1 ^= (s1 << 7) & 0x9d2c5680;
428  s1 ^= (s1 << 15) & 0xefc60000;
429  return (s1 ^ (s1 >> 18));
430 }
431 
432 inline double
434 {
435  return double(GetIntegerVariate()) * (1.0 / 4294967295.0);
436 }
437 
439 inline double
441 {
442  return GetVariateWithClosedRange() * n;
443 }
444 
446 inline double
448 {
449  return double(GetIntegerVariate()) * (1.0 / 4294967296.0);
450 }
451 
453 inline double
455 {
456  return GetVariateWithOpenUpperRange() * n;
457 }
458 
460 inline double
462 {
463  return (double(GetIntegerVariate()) + 0.5) * (1.0 / 4294967296.0);
464 }
465 
467 inline double
469 {
470  return GetVariateWithOpenRange() * n;
471 }
472 
475 {
476  // Find which bits are used in n
477  IntegerType used = n;
478 
479  used |= used >> 1;
480  used |= used >> 2;
481  used |= used >> 4;
482  used |= used >> 8;
483  used |= used >> 16;
484 
485  // Draw numbers until one is found in [0,n]
486  IntegerType i;
487  do
488  {
489  i = GetIntegerVariate() & used; // toss unused bits to shorten search
490  } while (i > n);
491 
492  return i;
493 }
494 
497 inline double
499 {
500  IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
501 
502  return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku
503  // Wada
504 }
505 
507 // TODO: Compare with vnl_sample_normal
508 inline double
509 MersenneTwisterRandomVariateGenerator::GetNormalVariate(const double & mean, const double & variance)
510 {
511  // Return a real number from a normal (Gaussian) distribution with given
512  // mean and variance by Box-Muller method
513  double r = std::sqrt(-2.0 * std::log(1.0 - GetVariateWithOpenRange()) * variance);
514  double phi = 2.0 * itk::Math::pi * GetVariateWithOpenUpperRange();
516 
517  return mean + r * std::cos(phi);
518 }
519 
521 // TODO: Compare with vnl_sample_uniform
522 inline double
524 {
525  double u = GetVariateWithOpenUpperRange();
526 
527  return ((1.0 - u) * a + u * b);
528 }
529 
530 inline double
532 {
533  return GetVariateWithClosedRange();
534 }
535 
536 inline double
538 {
539  return GetVariate();
540 }
541 
542 /* Change log from MTRand.h */
543 // Change log:
544 //
545 // v0.1 - First release on 15 May 2000
546 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
547 // - Translated from C to C++
548 // - Made completely ANSI compliant
549 // - Designed convenient interface for initialization, seeding, and
550 // obtaining numbers in default or user-defined ranges
551 // - Added automatic seeding from /dev/urandom or time() and clock()
552 // - Provided functions for saving and loading generator state
553 //
554 // v0.2 - Fixed bug which reloaded generator one step too late
555 //
556 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
557 //
558 // v0.4 - Removed trailing newline in saved generator format to be consistent
559 // with output format of built-in types
560 //
561 // v0.5 - Improved portability by replacing static const int's with enum's and
562 // clarifying return values in seed(); suggested by Eric Heimburg
563 // - Removed MAXINT constant; use 0xffffffffUL instead
564 //
565 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
566 // - Changed integer [0,n] generator to give better uniformity
567 //
568 // v0.7 - Fixed operator precedence ambiguity in reload()
569 // - Added access for real numbers in (0,1) and (0,n)
570 //
571 // v0.8 - Included time.h header to properly support time_t and clock_t
572 //
573 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
574 // - Allowed for seeding with arrays of any length
575 // - Added access for real numbers in [0,1) with 53-bit resolution
576 // - Added access for real numbers from normal (Gaussian) distributions
577 // - Increased overall speed by optimizing twist()
578 // - Doubled speed of integer [0,n] generation
579 // - Fixed out-of-range number generation on 64-bit machines
580 // - Improved portability by substituting literal constants for long enum's
581 // - Changed license from GNU LGPL to BSD
582 } // end namespace Statistics
583 } // end namespace itk
584 
585 #endif
itk::Statistics::MersenneTwisterRandomVariateGenerator::mixBits
IntegerType mixBits(const IntegerType &u, const IntegerType &v) const
Definition: itkMersenneTwisterRandomVariateGenerator.h:293
itk::Statistics::MersenneTwisterRandomVariateGenerator::hiBit
IntegerType hiBit(const IntegerType &u) const
Definition: itkMersenneTwisterRandomVariateGenerator.h:278
itkObjectFactory.h
itk::Statistics::MersenneTwisterRandomVariateGenerator
MersenneTwisterRandom random variate generator.
Definition: itkMersenneTwisterRandomVariateGenerator.h:127
itk::Statistics::MersenneTwisterRandomVariateGenerator::m_Left
int m_Left
Definition: itkMersenneTwisterRandomVariateGenerator.h:314
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetVariate
double GetVariate() override
Definition: itkMersenneTwisterRandomVariateGenerator.h:531
itk::Statistics::MersenneTwisterRandomVariateGenerator::hash
static IntegerType hash(time_t t, clock_t c)
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetNormalVariate
double GetNormalVariate(const double &mean=0.0, const double &variance=1.0)
Definition: itkMersenneTwisterRandomVariateGenerator.h:509
itk::Statistics::MersenneTwisterRandomVariateGenerator::state
IntegerType state[StateVectorLength]
Definition: itkMersenneTwisterRandomVariateGenerator.h:308
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetSeed
IntegerType GetSeed() const
Definition: itkMersenneTwisterRandomVariateGenerator.h:409
itkSingletonMacro.h
itk::Statistics::MersenneTwisterRandomVariateGenerator::loBit
IntegerType loBit(const IntegerType &u) const
Definition: itkMersenneTwisterRandomVariateGenerator.h:283
itk::Statistics::MersenneTwisterRandomVariateGenerator::SetSeed
void SetSeed()
Definition: itkMersenneTwisterRandomVariateGenerator.h:401
itk::SmartPointer< Self >
itk::Indent
Control indentation during Print() invocation.
Definition: itkIndent.h:49
itk::Statistics::MersenneTwisterRandomVariateGenerator::M
static constexpr unsigned int M
Definition: itkMersenneTwisterRandomVariateGenerator.h:271
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetUniformVariate
double GetUniformVariate(const double &a, const double &b)
Definition: itkMersenneTwisterRandomVariateGenerator.h:523
itk::Statistics::RandomVariateGeneratorBase
Defines common interfaces for random variate generators.
Definition: itkRandomVariateGeneratorBase.h:33
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange
double GetVariateWithOpenRange()
Definition: itkMersenneTwisterRandomVariateGenerator.h:461
itk::LightObject
Light weight base class for most itk classes.
Definition: itkLightObject.h:59
itkMacro.h
itkRandomVariateGeneratorBase.h
itk::Statistics::MersenneTwisterRandomVariateGenerator::m_InstanceLock
std::mutex m_InstanceLock
Definition: itkMersenneTwisterRandomVariateGenerator.h:328
itk::Statistics::MersenneTwisterRandomVariateGenerator::twist
IntegerType twist(const IntegerType &m, const IntegerType &s0, const IntegerType &s1) const
Definition: itkMersenneTwisterRandomVariateGenerator.h:299
itk::Statistics::MersenneTwisterRandomVariateGenerator::operator()
double operator()()
Definition: itkMersenneTwisterRandomVariateGenerator.h:537
itkIntTypes.h
itk::Statistics::MersenneTwisterRandomVariateGenerator::Get53BitVariate
double Get53BitVariate()
Definition: itkMersenneTwisterRandomVariateGenerator.h:498
itk::Statistics::MersenneTwisterRandomVariateGenerator::m_PNext
IntegerType * m_PNext
Definition: itkMersenneTwisterRandomVariateGenerator.h:311
itk::Statistics::MersenneTwisterRandomVariateGenerator::IntegerType
uint32_t IntegerType
Definition: itkMersenneTwisterRandomVariateGenerator.h:136
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange
double GetVariateWithClosedRange()
Definition: itkMersenneTwisterRandomVariateGenerator.h:433
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkArray.h:26
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange
double GetVariateWithOpenUpperRange()
Definition: itkMersenneTwisterRandomVariateGenerator.h:447
itkGetGlobalDeclarationMacro
#define itkGetGlobalDeclarationMacro(Type, VarName)
Definition: itkSingletonMacro.h:34
itk::Statistics::MersenneTwisterRandomVariateGenerator::m_Seed
std::atomic< IntegerType > m_Seed
Definition: itkMersenneTwisterRandomVariateGenerator.h:317
itk::Statistics::MersenneTwisterRandomVariateGenerator::loBits
IntegerType loBits(const IntegerType &u) const
Definition: itkMersenneTwisterRandomVariateGenerator.h:288
itk::uint32_t
::uint32_t uint32_t
Definition: itkIntTypes.h:33
itk::Statistics::MersenneTwisterRandomVariateGenerator::StateVectorLength
static constexpr IntegerType StateVectorLength
Definition: itkMersenneTwisterRandomVariateGenerator.h:164
itk::Math::pi
static constexpr double pi
Definition: itkMath.h:64
itkMath.h
itk::Statistics::MersenneTwisterRandomVariateGenerator::Initialize
void Initialize()
Definition: itkMersenneTwisterRandomVariateGenerator.h:388
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetIntegerVariate
IntegerType GetIntegerVariate()
Definition: itkMersenneTwisterRandomVariateGenerator.h:416
itk::Statistics::MersenneTwisterRandomVariateGenerator::m_PimplGlobals
static MersenneTwisterGlobals * m_PimplGlobals
Definition: itkMersenneTwisterRandomVariateGenerator.h:332
itk::Statistics::MersenneTwisterRandomVariateGenerator::reload
void reload()
Definition: itkMersenneTwisterRandomVariateGenerator.h:361