ITK  4.4.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 "vcl_ctime.h"
26 #include "vnl/vnl_math.h"
27 #include <climits>
28 
29 namespace itk
30 {
31 namespace Statistics
32 {
107 {
108 public:
109 
115 
117 
121 
125  static Pointer New();
126 
132  static Pointer GetInstance();
133 
135  itkStaticConstMacro(StateVectorLength, IntegerType, 624);
136 
138  void Initialize(const IntegerType oneSeed);
139 
140  /* Initialize with vcl_clock time */
141  void Initialize();
142 
144  double GetVariateWithClosedRange();
145 
147  double GetVariateWithClosedRange(const double & n);
148 
150  double GetVariateWithOpenUpperRange();
151 
153  double GetVariateWithOpenUpperRange(const double & n);
154 
156  double GetVariateWithOpenRange();
157 
159  double GetVariateWithOpenRange(const double & n);
160 
162  IntegerType GetIntegerVariate();
163 
165  IntegerType GetIntegerVariate(const IntegerType & n);
166 
169  double Get53BitVariate();
170 
171  /* Access to a normal random number distribution
172  * TODO: Compare with vnl_sample_normal */
173  double GetNormalVariate(const double & mean = 0.0,
174  const double & variance = 1.0);
175 
176  /* Access to a uniform random number distribution in the range [a, b)
177  * TODO: Compare with vnl_sample_uniform */
178  double GetUniformVariate(const double & a, const double & b);
179 
185  virtual double GetVariate();
186 
188  double operator()();
189 
190  // Re-seeding functions with same behavior as initializers
191  inline void SetSeed(const IntegerType oneSeed);
192  // inline void SetSeed(IntegerType *bigSeed, const IntegerType seedLength = StateVectorLength);
193  inline void SetSeed();
194  // Return the current seed
195  IntegerType GetSeed() { return this->m_Seed; };
196 
197  /*
198  // Saving and loading generator state
199  void save( IntegerType* saveArray ) const; // to array of size SAVE
200  void load( IntegerType *const loadArray ); // from such array
201  */
202 
203 protected:
206  virtual void PrintSelf(std::ostream & os, Indent indent) const
207  {
208  Superclass::PrintSelf(os, indent);
209 
210  // Print state vector contents
211  os << indent << "State vector: " << state << std::endl;
212  os << indent;
213  register const IntegerType *s = state;
214  register int i = StateVectorLength;
215  for (; i--; os << *s++ << "\t" ) {}
216  os << std::endl;
217 
218  //Print next value to be gotten from state
219  os << indent << "Next value to be gotten from state: " << pNext << std::endl;
220 
221  //Number of values left before reload
222  os << indent << "Values left before next reload: " << left << std::endl;
223  }
224 
225  // period parameter
226  itkStaticConstMacro (M, unsigned int, 397);
227 
228  IntegerType state[StateVectorLength]; // internal state
229  IntegerType *pNext; // next value to get from state
230  int left; // number of values left before reload
231  // needed
232  IntegerType m_Seed; // Seed value
233 
234  /* Reload array with N new values */
235  void reload();
236 
237  IntegerType hiBit(const IntegerType & u) const { return u & 0x80000000; }
238  IntegerType loBit(const IntegerType & u) const { return u & 0x00000001; }
239  IntegerType loBits(const IntegerType & u) const { return u & 0x7fffffff; }
240  IntegerType mixBits(const IntegerType & u, const IntegerType & v) const
241  {
242  return hiBit(u) | loBits(v);
243  }
244 
245  IntegerType twist(const IntegerType & m, const IntegerType & s0, const IntegerType & s1) const
246  {
247  return m ^ ( mixBits(s0, s1) >> 1 ) ^ ( -loBit(s1) & 0x9908b0df );
248  }
249 
250  static IntegerType hash(vcl_time_t t, vcl_clock_t c);
251 
253 }; // end of class
254 
255 // Declare inlined functions.... (must be declared in the header)
256 
259 {
260  // Get a IntegerType from t and c
261  // Better than IntegerType(x) in case x is floating point in [0,1]
262  // Based on code by Lawrence Kirby: fred at genesis dot demon dot co dot uk
263 
264  static IntegerType differ = 0; // guarantee time-based seeds will change
265 
266  IntegerType h1 = 0;
267  unsigned char *p = (unsigned char *)&t;
268 
269  const unsigned int sizeOfT = static_cast< unsigned int >( sizeof(t) );
270  for ( unsigned int i = 0; i < sizeOfT; ++i )
271  {
272  h1 *= UCHAR_MAX + 2U;
273  h1 += p[i];
274  }
275  IntegerType h2 = 0;
276  p = (unsigned char *)&c;
277 
278  const unsigned int sizeOfC = static_cast< unsigned int >( sizeof(c) );
279  for ( unsigned int j = 0; j < sizeOfC; ++j )
280  {
281  h2 *= UCHAR_MAX + 2U;
282  h2 += p[j];
283  }
284  return ( h1 + differ++ ) ^ h2;
285 }
286 
287 inline void
289 {
290  this->m_Seed = seed;
291  // Initialize generator state with seed
292  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
293  // In previous versions, most significant bits (MSBs) of the seed affect
294  // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
295  register IntegerType *s = state;
296  register IntegerType *r = state;
297  register IntegerType i = 1;
298 
299  *s++ = seed & 0xffffffffUL;
301  {
302  *s++ = ( 1812433253UL * ( *r ^ ( *r >> 30 ) ) + i ) & 0xffffffffUL;
303  r++;
304  }
305 }
306 
307 inline void
309 {
310  // Generate N new values in state
311  // Made clearer and faster by Matthew Bellew
312  // matthew dot bellew at home dot com
313 
314  // get rid of VS warning
315  register int index = static_cast< int >(
317 
318  register IntegerType *p = state;
319  register int i;
320 
322  {
323  *p = twist(p[M], p[0], p[1]);
324  }
325  for ( i = M; --i; ++p )
326  {
327  *p = twist(p[index], p[0], p[1]);
328  }
329  *p = twist(p[index], p[0], state[0]);
330 
332 }
333 
334  /*
335 #define SVL 624
336 inline void
337 MersenneTwisterRandomVariateGenerator::SetSeed(
338  IntegerType *const bigSeed, const IntegerType seedLength)
339 {
340  // Seed the generator with an array of IntegerType's
341  // There are 2^19937-1 possible initial states. This function allows
342  // all of those to be accessed by providing at least 19937 bits (with a
343  // default seed length of StateVectorLength = 624 IntegerType's).
344  // Any bits above the lower 32
345  // in each element are discarded.
346  // Just call seed() if you want to get array from /dev/urandom
347  Initialize(19650218UL);
348  register IntegerType i = 1;
349  register IntegerType j = 0;
350  register int k;
351  if ( StateVectorLength > seedLength )
352  {
353  k = StateVectorLength;
354  }
355  else
356  {
357  k = seedLength;
358  }
359  for (; k; --k )
360  {
361  state[i] =
362  state[i] ^ ( ( state[i - 1] ^ ( state[i - 1] >> 30 ) ) * 1664525UL );
363  state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
364  state[i] &= 0xffffffffUL;
365  ++i; ++j;
366  if ( i >= StateVectorLength ) { state[0] = state[StateVectorLength - 1]; i = 1; }
367  if ( j >= seedLength ) { j = 0; }
368  }
369  for ( k = StateVectorLength - 1; k; --k )
370  {
371  state[i] =
372  state[i] ^ ( ( state[i - 1] ^ ( state[i - 1] >> 30 ) ) * 1566083941UL );
373  state[i] -= i;
374  state[i] &= 0xffffffffUL;
375  ++i;
376  if ( i >= SVL )
377  {
378  state[0] = state[StateVectorLength - 1]; i = 1;
379  }
380  }
381  state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
382  reload();
383 }
384  */
385 
386 inline void
388 {
389  SetSeed();
390 }
391 
392 inline void
394 {
395  // Seed the generator with a simple IntegerType
396  Initialize(oneSeed);
397  reload();
398 }
399 
400 inline void
402 {
403  // use time() and clock() to generate a unlikely-to-repeat seed.
404  SetSeed( hash( vcl_time(0), vcl_clock() ) );
405 }
406 
410 {
411  if ( left == 0 ) { reload(); }
412  --left;
414 
415  register IntegerType s1;
416  s1 = *pNext++;
417  s1 ^= ( s1 >> 11 );
418  s1 ^= ( s1 << 7 ) & 0x9d2c5680;
419  s1 ^= ( s1 << 15 ) & 0xefc60000;
420  return ( s1 ^ ( s1 >> 18 ) );
421 }
422 
423 inline double
425 {
426  return double( GetIntegerVariate() ) * ( 1.0 / 4294967295.0 );
427 }
428 
430 inline double
432  const double & n)
433 {
434  return GetVariateWithClosedRange() * n;
435 }
436 
438 inline double
440 {
441  return double( GetIntegerVariate() ) * ( 1.0 / 4294967296.0 );
442 }
443 
445 inline double
447  const double & n)
448 {
449  return GetVariateWithOpenUpperRange() * n;
450 }
451 
453 inline double
455 {
456  return ( double( GetIntegerVariate() ) + 0.5 ) * ( 1.0 / 4294967296.0 );
457 }
458 
460 inline double
462  const double & n)
463 {
464  return GetVariateWithOpenRange() * n;
465 }
466 
469  const IntegerType & n)
470 {
471  // Find which bits are used in n
472  // Optimized by Magnus Jonsson magnus at smartelectronix dot com
473  IntegerType used = n;
474 
475  used |= used >> 1;
476  used |= used >> 2;
477  used |= used >> 4;
478  used |= used >> 8;
479  used |= used >> 16;
480 
481  // Draw numbers until one is found in [0,n]
482  IntegerType i;
483  do
484  {
485  i = GetIntegerVariate() & used; // toss unused bits to shorten search
486  }
487  while ( i > n );
488 
489  return i;
490 }
491 
494 inline double
496 {
497  IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
498 
499  return ( a * 67108864.0 + b ) * ( 1.0 / 9007199254740992.0 ); // by Isaku
500  // Wada
501 }
502 
503 /* Access to a normal random number distribution */
504 // TODO: Compare with vnl_sample_normal
505 inline double
507  const double & mean, const double & variance)
508 {
509  // Return a real number from a normal (Gaussian) distribution with given
510  // mean and variance by Box-Muller method
511  double r = vcl_sqrt(-2.0 * vcl_log( 1.0 - GetVariateWithOpenRange() ) * variance);
512  double phi = 2.0 * vnl_math::pi
514 
515  return mean + r *vcl_cos(phi);
516 }
517 
518 /* Access to a uniform random number distribution */
519 // TODO: Compare with vnl_sample_uniform
520 inline double
522  const double & a, const double & b)
523 {
524  double u = GetVariateWithOpenUpperRange();
525 
526  return ( ( 1.0 - u ) * a + u * b );
527 }
528 
529 inline double
531 {
532  return GetVariateWithClosedRange();
533 }
534 
535 inline double
537 {
538  return GetVariate();
539 }
540 
541 inline
543 {
544  SetSeed (121212);
545 }
546 
547 /* Change log from MTRand.h */
548 // Change log:
549 //
550 // v0.1 - First release on 15 May 2000
551 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
552 // - Translated from C to C++
553 // - Made completely ANSI compliant
554 // - Designed convenient interface for initialization, seeding, and
555 // obtaining numbers in default or user-defined ranges
556 // - Added automatic seeding from /dev/urandom or time() and clock()
557 // - Provided functions for saving and loading generator state
558 //
559 // v0.2 - Fixed bug which reloaded generator one step too late
560 //
561 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
562 //
563 // v0.4 - Removed trailing newline in saved generator format to be consistent
564 // with output format of built-in types
565 //
566 // v0.5 - Improved portability by replacing static const int's with enum's and
567 // clarifying return values in seed(); suggested by Eric Heimburg
568 // - Removed MAXINT constant; use 0xffffffffUL instead
569 //
570 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
571 // - Changed integer [0,n] generator to give better uniformity
572 //
573 // v0.7 - Fixed operator precedence ambiguity in reload()
574 // - Added access for real numbers in (0,1) and (0,n)
575 //
576 // v0.8 - Included time.h header to properly support time_t and clock_t
577 //
578 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
579 // - Allowed for seeding with arrays of any length
580 // - Added access for real numbers in [0,1) with 53-bit resolution
581 // - Added access for real numbers from normal (Gaussian) distributions
582 // - Increased overall speed by optimizing twist()
583 // - Doubled speed of integer [0,n] generation
584 // - Fixed out-of-range number generation on 64-bit machines
585 // - Improved portability by substituting literal constants for long enum's
586 // - Changed license from GNU LGPL to BSD
587 } // end of namespace Statistics
588 } // end of namespace itk
589 
590 #endif
591