ITK  4.9.0
Insight Segmentation and Registration 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 "vnl/vnl_math.h"
26 #include <climits>
27 #include <ctime>
28 
29 namespace itk
30 {
31 namespace Statistics
32 {
109 {
110 public:
111 
117 
119 
123 
129  static Pointer New();
130 
137  static Pointer GetInstance();
138 
140  itkStaticConstMacro(StateVectorLength, IntegerType, 624);
141 
143  void Initialize(const IntegerType oneSeed);
144 
145  /* Initialize with clock time */
146  void Initialize();
147 
149  double GetVariateWithClosedRange();
150 
152  double GetVariateWithClosedRange(const double & n);
153 
155  double GetVariateWithOpenUpperRange();
156 
158  double GetVariateWithOpenUpperRange(const double & n);
159 
161  double GetVariateWithOpenRange();
162 
164  double GetVariateWithOpenRange(const double & n);
165 
167  IntegerType GetIntegerVariate();
168 
170  IntegerType GetIntegerVariate(const IntegerType & n);
171 
174  double Get53BitVariate();
175 
176  /* Access to a normal random number distribution
177  * TODO: Compare with vnl_sample_normal */
178  double GetNormalVariate(const double & mean = 0.0,
179  const double & variance = 1.0);
180 
181  /* Access to a uniform random number distribution in the range [a, b)
182  * TODO: Compare with vnl_sample_uniform */
183  double GetUniformVariate(const double & a, const double & b);
184 
190  virtual double GetVariate() ITK_OVERRIDE;
191 
193  double operator()();
194 
196  inline void SetSeed(const IntegerType oneSeed);
197  inline void SetSeed();
199 
201  IntegerType GetSeed() { return this->m_Seed; };
202 
203  /*
204  // Saving and loading generator state
205  void save( IntegerType* saveArray ) const; // to array of size SAVE
206  void load( IntegerType *const loadArray ); // from such array
207  */
208 
209 protected:
212  virtual void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE;
213 
214  // period parameter
215  itkStaticConstMacro (M, unsigned int, 397);
216 
217  IntegerType state[StateVectorLength]; // internal state
218  IntegerType *pNext; // next value to get from state
219  int left; // number of values left before reload
220  // needed
221  IntegerType m_Seed; // Seed value
222 
223  /* Reload array with N new values */
224  void reload();
225 
226  IntegerType hiBit(const IntegerType & u) const { return u & 0x80000000; }
227  IntegerType loBit(const IntegerType & u) const { return u & 0x00000001; }
228  IntegerType loBits(const IntegerType & u) const { return u & 0x7fffffff; }
229  IntegerType mixBits(const IntegerType & u, const IntegerType & v) const
230  {
231  return hiBit(u) | loBits(v);
232  }
233 
234  IntegerType twist(const IntegerType & m, const IntegerType & s0, const IntegerType & s1) const
235  {
236  return m ^ ( mixBits(s0, s1) >> 1 ) ^ ( -static_cast<int32_t>(loBit(s1)) & 0x9908b0df );
237  }
238 
239  static IntegerType hash(time_t t, clock_t c);
240 
241 private:
242 
244  static Pointer CreateInstance();
245 
246  // Static/Global Variable need to be thread-safely accessed
250 
251 }; // end of class
252 
253 // Declare inlined functions.... (must be declared in the header)
254 
255 inline void
257 {
258  this->m_Seed = seed;
259  // Initialize generator state with seed
260  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
261  // In previous versions, most significant bits (MSBs) of the seed affect
262  // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
263  IntegerType *s = state;
264  IntegerType *r = state;
265  IntegerType i = 1;
266 
267  *s++ = seed & 0xffffffffUL;
269  {
270  *s++ = ( 1812433253UL * ( *r ^ ( *r >> 30 ) ) + i ) & 0xffffffffUL;
271  r++;
272  }
273 }
274 
275 inline void
277 {
278  // Generate N new values in state
279  // Made clearer and faster by Matthew Bellew
280  // matthew dot bellew at home dot com
281 
282  // get rid of VS warning
283  int index = static_cast< int >(
285 
286  IntegerType *p = state;
287  int i;
288 
290  {
291  *p = twist(p[M], p[0], p[1]);
292  }
293  for ( i = M; --i; ++p )
294  {
295  *p = twist(p[index], p[0], p[1]);
296  }
297  *p = twist(p[index], p[0], state[0]);
298 
300 }
301 
302 inline void
304 {
305  SetSeed();
306 }
307 
308 inline void
310 {
311  // Seed the generator with a simple IntegerType
312  Initialize(oneSeed);
313  reload();
314 }
315 
316 inline void
318 {
319  // use time() and clock() to generate a unlikely-to-repeat seed.
320  SetSeed( hash( time(ITK_NULLPTR), clock() ) );
321 }
322 
326 {
327  if ( left == 0 ) { reload(); }
328  --left;
330 
331  IntegerType s1 = *pNext++;
332  s1 ^= ( s1 >> 11 );
333  s1 ^= ( s1 << 7 ) & 0x9d2c5680;
334  s1 ^= ( s1 << 15 ) & 0xefc60000;
335  return ( s1 ^ ( s1 >> 18 ) );
336 }
337 
338 inline double
340 {
341  return double( GetIntegerVariate() ) * ( 1.0 / 4294967295.0 );
342 }
343 
345 inline double
347  const double & n)
348 {
349  return GetVariateWithClosedRange() * n;
350 }
351 
353 inline double
355 {
356  return double( GetIntegerVariate() ) * ( 1.0 / 4294967296.0 );
357 }
358 
360 inline double
362  const double & n)
363 {
364  return GetVariateWithOpenUpperRange() * n;
365 }
366 
368 inline double
370 {
371  return ( double( GetIntegerVariate() ) + 0.5 ) * ( 1.0 / 4294967296.0 );
372 }
373 
375 inline double
377  const double & n)
378 {
379  return GetVariateWithOpenRange() * n;
380 }
381 
384  const IntegerType & n)
385 {
386  // Find which bits are used in n
387  // Optimized by Magnus Jonsson magnus at smartelectronix dot com
388  IntegerType used = n;
389 
390  used |= used >> 1;
391  used |= used >> 2;
392  used |= used >> 4;
393  used |= used >> 8;
394  used |= used >> 16;
395 
396  // Draw numbers until one is found in [0,n]
397  IntegerType i;
398  do
399  {
400  i = GetIntegerVariate() & used; // toss unused bits to shorten search
401  }
402  while ( i > n );
403 
404  return i;
405 }
406 
409 inline double
411 {
412  IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
413 
414  return ( a * 67108864.0 + b ) * ( 1.0 / 9007199254740992.0 ); // by Isaku
415  // Wada
416 }
417 
418 /* Access to a normal random number distribution */
419 // TODO: Compare with vnl_sample_normal
420 inline double
422  const double & mean, const double & variance)
423 {
424  // Return a real number from a normal (Gaussian) distribution with given
425  // mean and variance by Box-Muller method
426  double r = std::sqrt(-2.0 * std::log( 1.0 - GetVariateWithOpenRange() ) * variance);
427  double phi = 2.0 * vnl_math::pi
429 
430  return mean + r *std::cos(phi);
431 }
432 
433 /* Access to a uniform random number distribution */
434 // TODO: Compare with vnl_sample_uniform
435 inline double
437  const double & a, const double & b)
438 {
439  double u = GetVariateWithOpenUpperRange();
440 
441  return ( ( 1.0 - u ) * a + u * b );
442 }
443 
444 inline double
446 {
447  return GetVariateWithClosedRange();
448 }
449 
450 inline double
452 {
453  return GetVariate();
454 }
455 
456 /* Change log from MTRand.h */
457 // Change log:
458 //
459 // v0.1 - First release on 15 May 2000
460 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
461 // - Translated from C to C++
462 // - Made completely ANSI compliant
463 // - Designed convenient interface for initialization, seeding, and
464 // obtaining numbers in default or user-defined ranges
465 // - Added automatic seeding from /dev/urandom or time() and clock()
466 // - Provided functions for saving and loading generator state
467 //
468 // v0.2 - Fixed bug which reloaded generator one step too late
469 //
470 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
471 //
472 // v0.4 - Removed trailing newline in saved generator format to be consistent
473 // with output format of built-in types
474 //
475 // v0.5 - Improved portability by replacing static const int's with enum's and
476 // clarifying return values in seed(); suggested by Eric Heimburg
477 // - Removed MAXINT constant; use 0xffffffffUL instead
478 //
479 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
480 // - Changed integer [0,n] generator to give better uniformity
481 //
482 // v0.7 - Fixed operator precedence ambiguity in reload()
483 // - Added access for real numbers in (0,1) and (0,n)
484 //
485 // v0.8 - Included time.h header to properly support time_t and clock_t
486 //
487 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
488 // - Allowed for seeding with arrays of any length
489 // - Added access for real numbers in [0,1) with 53-bit resolution
490 // - Added access for real numbers from normal (Gaussian) distributions
491 // - Increased overall speed by optimizing twist()
492 // - Doubled speed of integer [0,n] generation
493 // - Fixed out-of-range number generation on 64-bit machines
494 // - Improved portability by substituting literal constants for long enum's
495 // - Changed license from GNU LGPL to BSD
496 } // end of namespace Statistics
497 } // end of namespace itk
498 
499 #endif
Critical section locking class that can be allocated on the stack.
static const double pi
Definition: itkMath.h:57
KWIML_INT_uint32_t uint32_t
Definition: itkIntTypes.h:87
IntegerType mixBits(const IntegerType &u, const IntegerType &v) const
Defines common interfaces for random variate generators.
IntegerType twist(const IntegerType &m, const IntegerType &s0, const IntegerType &s1) const
double GetNormalVariate(const double &mean=0.0, const double &variance=1.0)
Control indentation during Print() invocation.
Definition: itkIndent.h:49
static IntegerType hash(time_t t, clock_t c)