KWStyle - itkMersenneTwisterRandomVariateGenerator.h

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkMersenneTwisterRandomVariateGenerator.h.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:41 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17		#ifndef __itkMersenneTwisterRandomVariateGenerator_h
18		#define __itkMersenneTwisterRandomVariateGenerator_h
19
20		#include "itkMacro.h"
21		#include "itkObjectFactory.h"
22		#include "itkRandomVariateGeneratorBase.h"
23		#include "itkIntTypes.h"
24		#include "vcl_ctime.h"
25		#include "vnl/vnl_math.h"
26		#include <limits.h>
27
28
29		namespace itk {
30		namespace Statistics {
31
32		/** \class MersenneTwisterRandomVariateGenerator
33		* \brief MersenneTwisterRandom random variate generator
34		*
35		* This notice was included with the original implementation.
36		* The only changes made were to obfuscate the author's email addresses.
37		*
38		* MersenneTwister.h
39		* Mersenne Twister random number generator -- a C++ class MTRand
40		* Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
41		* Richard J. Wagner  v1.0  15 May 2003  rjwagner at writeme dot com
42		*
43		* The Mersenne Twister is an algorithm for generating random numbers.  It
44		* was designed with consideration of the flaws in various other generators.
45		* The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
46		* are far greater.  The generator is also fast; it avoids multiplication and
47		* division, and it benefits from caches and pipelines.  For more information
48		* see the inventors' web page at http:*www.math.keio.ac.jp/~matumoto/emt.html
49
50		* Reference
51		* M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
52		* Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
53		* Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
54		*
55		* Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
56		* Copyright (C) 2000 - 2003, Richard J. Wagner
57		* All rights reserved.
58		*
59		* Redistribution and use in source and binary forms, with or without
60		* modification, are permitted provided that the following conditions
61		* are met:
62		*
63		*   1. Redistributions of source code must retain the above copyright
64		*      notice, this list of conditions and the following disclaimer.
65		*
66		*   2. Redistributions in binary form must reproduce the above copyright
67		*      notice, this list of conditions and the following disclaimer in the
68		*      documentation and/or other materials provided with the distribution.
69		*
70		*   3. The names of its contributors may not be used to endorse or promote
71		*      products derived from this software without specific prior written
72		*      permission.
73		*
74		* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
75		* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
76		* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
77	LEN	* A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
78		* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
79		* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
80		* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
81		* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
82		* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
83		* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
84		* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
85		*
86		* The original code included the following notice:
87		*
88	LEN	*     When you use this, send an email to: matumoto at math dot keio dot ac dot jp
89		*     with an appropriate reference to your work.
90		*
91		* It would be nice to CC:
92		* rjwagner at writeme dot com and Cokus at math dot washington dot edu
93		* when you write.
94		*
95		* \ingroup Statistics
96		*/
97		class ITKCommon_EXPORT MersenneTwisterRandomVariateGenerator :
98	IND	****public RandomVariateGeneratorBase
99		{
100		public:
101
102		/** Standard class typedefs. */
103	SEM	typedef MersenneTwisterRandomVariateGenerator Self ;
104	TDA	typedef RandomVariateGeneratorBase Superclass;
105	TDA	typedef SmartPointer<Self>   Pointer;
106	TDA	typedef SmartPointer<const Self>  ConstPointer;
107	TDA	typedef ITK_UINT32 IntegerType;
108
109		/** Run-time type information (and related methods). */
110		itkTypeMacro(MersenneTwisterRandomVariateGenerator,
111		RandomVariateGeneratorBase );
112
113		/** Method for creation through the object factory. */
114		static Pointer New();
115		static Pointer GetInstance();
116
117		/** Length of state vector */
118		itkStaticConstMacro(StateVectorLength,IntegerType,624);
119
120		/** Length of array for save() */
121		// itkStaticConstMacro(SAVE,IntegerType,625);
122
123
124
125		// Methods to reseed
126
127		/** initialize with a simple IntegerType */
128		void Initialize( const IntegerType oneSeed );
129
130		/** Initialize with an array */
131		//void Initialize( IntegerType *const bigSeed,
132		//        IntegerType const seedLength = StateVectorLength );
133
134		/* Initialize with vcl_clock time */
135		void Initialize();
136
137		// Methods to get various random variates ...
138
139		/** Get a random variate in the range [0, 1] */
140		double GetVariateWithClosedRange();
141
142		/** Get a random variate in the range [0, n] */
143		double GetVariateWithClosedRange( const double& n );
144
145		/** Get a range variate in the range [0, 1) */
146		double GetVariateWithOpenUpperRange();
147
148		/** Get a range variate in the range [0, n) */
149		double GetVariateWithOpenUpperRange( const double& n );
150
151		/** Get a range variate in the range (0, 1) */
152		double GetVariateWithOpenRange();
153
154		/** Get a range variate in the range (0, n) */
155		double GetVariateWithOpenRange( const double& n );
156
157		/** Get an integer variate in [0, 2^32-1] */
158		IntegerType GetIntegerVariate();
159
160		/** Get an integer variate in [0, n] for n < 2^32 */
161		IntegerType GetIntegerVariate( const IntegerType& n );
162
163		/** Access to 53-bit random numbers (capacity of IEEE double precision)
164		* in the range [0,1) */
165		double Get53BitVariate();
166
167		/* Access to a normal random number distribution
168		* TODO: Compare with vnl_sample_normal */
169		double GetNormalVariate( const double& mean = 0.0,
170		const double& variance = 1.0 );
171
172		/* Access to a uniform random number distribution in the range [a, b)
173		* TODO: Compare with vnl_sample_uniform */
174		double GetUniformVariate( const double& a, const double& b );
175
176		/** get a variate in the range [0, 1]
177		* Do NOT use for CRYPTOGRAPHY without securely hashing several returned
178		* values together, otherwise the generator state can be learned after
179		* reading 624 consecutive values.
180		*/
181		virtual double GetVariate();
182
183		/** Same as GetVariate() */
184		double operator()();
185
186
187		// Re-seeding functions with same behavior as initializers
188		inline void SetSeed( const IntegerType oneSeed );
189	LEN	inline void SetSeed( IntegerType * bigSeed, const IntegerType seedLength = StateVectorLength );
190		inline void SetSeed();
191
192		/*
193		// Saving and loading generator state
194		void save( IntegerType* saveArray ) const;  // to array of size SAVE
195		void load( IntegerType *const loadArray );  // from such array
196		*/
197
198	IND	*protected:
199		inline MersenneTwisterRandomVariateGenerator();
200		virtual ~MersenneTwisterRandomVariateGenerator() {};
201		virtual void PrintSelf(std::ostream& os, Indent indent) const {
202		Superclass::PrintSelf(os, indent);
203
204		// Print state vector contents
205		os << indent << "State vector: " << state << std::endl;
206		os << indent;
207		register const IntegerType *s = state;
208		register int i = StateVectorLength;
209	SEM	for( ; i--; os << *s++ << "\t" ) {}
210		os << std::endl;
211
212		//Print next value to be gotten from state
213	LEN	os << indent << "Next value to be gotten from state: " << pNext << std::endl;
214
215		//Number of values left before reload
216		os << indent << "Values left before next reload: " << left << std::endl;
217		}
218
219
220		// enum { M = 397 };  // period parameter
221		itkStaticConstMacro ( M, unsigned int, 397 );
222
223		IntegerType state[StateVectorLength];   // internal state
224		IntegerType *pNext;     // next value to get from state
225		int left;          // number of values left before reload needed
226
227		/* Reload array with N new values */
228		void reload();
229
230		IntegerType hiBit( const IntegerType& u ) const { return u & 0x80000000UL; }
231		IntegerType loBit( const IntegerType& u ) const { return u & 0x00000001UL; }
232		IntegerType loBits( const IntegerType& u ) const { return u & 0x7fffffffUL; }
233		IntegerType mixBits( const IntegerType& u, const IntegerType& v ) const
234		{
235		return hiBit(u) | loBits(v);
236		}
237	LEN	IntegerType twist( const IntegerType& m, const IntegerType& s0, const IntegerType& s1 ) const
238		{
239		return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL);
240		}
241		static IntegerType hash( vcl_time_t t, vcl_clock_t c );
242		static Pointer m_Instance;
243	SEM	} ;  // end of class
244
245
246
247	EML
248		// Declare inlined functions.... (must be declared in the header)
249		// Declare then in order to keep SGI happy
250
251		inline MersenneTwisterRandomVariateGenerator::IntegerType
252	IND	**MersenneTwisterRandomVariateGenerator::hash( vcl_time_t t, vcl_clock_t c )
253	IND	**{
254		// Get a IntegerType from t and c
255		// Better than IntegerType(x) in case x is floating point in [0,1]
256		// Based on code by Lawrence Kirby: fred at genesis dot demon dot co dot uk
257
258		static IntegerType differ = 0;  // guarantee time-based seeds will change
259
260		IntegerType h1 = 0;
261		unsigned char *p = (unsigned char *) &t;
262		for( size_t i = 0; i < sizeof(t); ++i )
263		{
264		h1 *= UCHAR_MAX + 2U;
265		h1 += p[i];
266		}
267		IntegerType h2 = 0;
268		p = (unsigned char *) &c;
269		for( size_t j = 0; j < sizeof(c); ++j )
270		{
271		h2 *= UCHAR_MAX + 2U;
272		h2 += p[j];
273		}
274		return ( h1 + differ++ ) ^ h2;
275	IND	**}
276
277
278		inline void
279		MersenneTwisterRandomVariateGenerator::Initialize( const IntegerType seed )
280	IND	**{
281		// Initialize generator state with seed
282		// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
283		// In previous versions, most significant bits (MSBs) of the seed affect
284		// only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
285		register IntegerType *s = state;
286		register IntegerType *r = state;
287		register IntegerType i = 1;
288		*s++ = seed & 0xffffffffUL;
289	LEN	for( i = 1; i < MersenneTwisterRandomVariateGenerator::StateVectorLength; ++i )
290		{
291		*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
292		r++;
293		}
294	IND	**}
295
296		inline void
297	IND	**MersenneTwisterRandomVariateGenerator::reload()
298	IND	**{
299		// Generate N new values in state
300		// Made clearer and faster by Matthew Bellew
301		// matthew dot bellew at home dot com
302
303		// get rid of VS warning
304		register int index = static_cast< int >(
305		M-MersenneTwisterRandomVariateGenerator::StateVectorLength);
306
307		register IntegerType *p = state;
308		register int i;
309	LEN	for( i = MersenneTwisterRandomVariateGenerator::StateVectorLength - M; i--; ++p )
310		{
311		*p = twist( p[M], p[0], p[1] );
312		}
313		for( i = M; --i; ++p )
314		{
315		*p = twist( p[index], p[0], p[1] );
316		}
317		*p = twist( p[index], p[0], state[0] );
318
319	LEN	left = MersenneTwisterRandomVariateGenerator::StateVectorLength, pNext = state;
320	IND	**}
321
322
323	EML
324		#define SVL 624
325		inline void
326	IND	**MersenneTwisterRandomVariateGenerator::SetSeed(
327		IntegerType * const bigSeed, const IntegerType seedLength )
328	IND	**{
329		// Seed the generator with an array of IntegerType's
330		// There are 2^19937-1 possible initial states.  This function allows
331		// all of those to be accessed by providing at least 19937 bits (with a
332		// default seed length of StateVectorLength = 624 IntegerType's).
333		// Any bits above the lower 32
334		// in each element are discarded.
335		// Just call seed() if you want to get array from /dev/urandom
336		Initialize(19650218UL);
337		register IntegerType i = 1;
338		register IntegerType j = 0;
339		register int k;
340		if ( StateVectorLength > seedLength )
341		{
342		k = StateVectorLength;
343		}
344		else
345		{
346		k = seedLength;
347		}
348	SEM	for( ; k; --k )
349		{
350		state[i] =
351	IND	******state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
352		state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
353		state[i] &= 0xffffffffUL;
354		++i;  ++j;
355	LEN	if( i >= StateVectorLength ) { state[0] = state[StateVectorLength-1];  i = 1; }
356		if( j >= seedLength ) j = 0;
357		}
358		for( k = StateVectorLength - 1; k; --k )
359		{
360		state[i] =
361	IND	******state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
362		state[i] -= i;
363		state[i] &= 0xffffffffUL;
364		++i;
365		if( i >= SVL )
366		{
367		state[0] = state[StateVectorLength-1];  i = 1;
368		}
369		}
370		state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
371		reload();
372		}
373
374
375		inline void
376	IND	**MersenneTwisterRandomVariateGenerator::Initialize()
377	IND	**{
378		SetSeed();
379	IND	**}
380
381
382		inline void
383	IND	**MersenneTwisterRandomVariateGenerator::SetSeed( const IntegerType oneSeed )
384	IND	**{
385		// Seed the generator with a simple IntegerType
386		Initialize(oneSeed);
387		reload();
388	IND	**}
389
390
391		inline void
392	IND	**MersenneTwisterRandomVariateGenerator::SetSeed()
393	IND	**{
394		/*
395		// Seed the generator with an array from /dev/urandom if available
396		// Otherwise use a hash of time() and clock() values
397
398		// First try getting an array from /dev/urandom
399		FILE* urandom = fopen( "/dev/urandom", "rb" );
400		if( urandom )
401		{
402	LEN	IntegerType bigSeed[MersenneTwisterRandomVariateGenerator::StateVectorLength];
403		register IntegerType *s = bigSeed;
404	LEN	register IntegerType i = MersenneTwisterRandomVariateGenerator::StateVectorLength;
405		register bool success = true;
406		while( success && i-- )
407	IND	******success = fread( s++, sizeof(IntegerType), 1, urandom );
408		fclose(urandom);
409	LEN	if( success ) { seed( bigSeed, MersenneTwisterRandomVariateGenerator::StateVectorLength );  return; }
410		}
411	IND	****/
412	IND	***// Was not successful, so use time() and clock() instead
413	IND	***// /dev/urandom is not portable, just use the clock
414		SetSeed( hash( vcl_time(0), vcl_clock() ) );
415	IND	**}
416
417
418	EML
419		/** Get an integer variate in [0, 2^32-1] */
420		inline MersenneTwisterRandomVariateGenerator::IntegerType
421	IND	**MersenneTwisterRandomVariateGenerator::GetIntegerVariate()
422	IND	**{
423		if( left == 0 ) reload();
424		--left;
425
426		register IntegerType s1;
427		s1 = *pNext++;
428		s1 ^= (s1 >> 11);
429		s1 ^= (s1 <<  7) & 0x9d2c5680UL;
430		s1 ^= (s1 << 15) & 0xefc60000UL;
431		return ( s1 ^ (s1 >> 18) );
432	IND	**}
433
434
435		inline double
436	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange()
437	IND	**{
438		return double(GetIntegerVariate()) * (1.0/4294967295.0);
439	IND	**}
440
441		/** Get a random variate in the range [0, n] */
442		inline double
443	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithClosedRange(
444		const double& n )
445	IND	**{
446		return rand() * n;
447	IND	**}
448
449		/** Get a range variate in the range [0, 1) */
450		inline double
451	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange()
452	IND	**{
453		return double(GetIntegerVariate()) * (1.0/4294967296.0);
454	IND	**}
455
456		/** Get a range variate in the range [0, n) */
457		inline double
458	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithOpenUpperRange(
459		const double& n )
460	IND	**{
461		return GetVariateWithOpenUpperRange() * n;
462	IND	**}
463
464		/** Get a range variate in the range (0, 1) */
465		inline double
466	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange()
467	IND	**{
468		return ( double(GetIntegerVariate()) + 0.5 ) * (1.0/4294967296.0);
469	IND	**}
470
471
472		/** Get a range variate in the range (0, n) */
473		inline double
474	IND	**MersenneTwisterRandomVariateGenerator::GetVariateWithOpenRange(
475		const double& n )
476	IND	**{
477		return GetVariateWithOpenRange() * n;
478	IND	**}
479
480
481		inline MersenneTwisterRandomVariateGenerator::IntegerType
482	IND	**MersenneTwisterRandomVariateGenerator::GetIntegerVariate(
483		const IntegerType& n )
484	IND	**{
485		// Find which bits are used in n
486		// Optimized by Magnus Jonsson magnus at smartelectronix dot com
487		IntegerType used = n;
488		used |= used >> 1;
489		used |= used >> 2;
490		used |= used >> 4;
491		used |= used >> 8;
492		used |= used >> 16;
493
494		// Draw numbers until one is found in [0,n]
495		IntegerType i;
496		do
497		{
498		i = GetIntegerVariate() & used;  // toss unused bits to shorten search
499		}  while( i > n );
500
501		return i;
502	IND	**}
503
504
505	EML
506		/** Access to 53-bit random numbers (capacity of IEEE double precision)
507		* in the range [0,1) */
508		inline double
509	IND	**MersenneTwisterRandomVariateGenerator::Get53BitVariate()
510	IND	**{
511		IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
512		return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
513	IND	**}
514
515
516		/* Access to a normal randon number distribution */
517		// TODO: Compare with vnl_sample_normal
518		inline double
519	IND	**MersenneTwisterRandomVariateGenerator::GetNormalVariate(
520		const double& mean, const double& variance )
521	IND	**{
522		// Return a real number from a normal (Gaussian) distribution with given
523		// mean and variance by Box-Muller method
524	LEN	double r = vcl_sqrt( -2.0 * vcl_log( 1.0-GetVariateWithOpenRange()) ) * variance;
525		double phi = 2.0 * 3.14159265358979323846264338328
526	IND	*************************** GetVariateWithOpenUpperRange();
527		return mean + r * vcl_cos(phi);
528	IND	**}
529
530
531	EML
532		/* Access to a uniform random number distribution */
533		// TODO: Compare with vnl_sample_uniform
534		inline double
535	IND	**MersenneTwisterRandomVariateGenerator::GetUniformVariate(
536		const double& a, const double& b )
537	IND	**{
538		double u = GetVariateWithOpenUpperRange();
539		return ((1.0 -u) * a + u * b);
540	IND	**}
541
542
543		inline double
544	IND	**MersenneTwisterRandomVariateGenerator::GetVariate()
545	IND	**{
546		return GetVariateWithClosedRange();
547	IND	**}
548
549
550		inline double
551	IND	**MersenneTwisterRandomVariateGenerator::operator()()
552	IND	**{
553		return GetVariate();
554	IND	**}
555
556
557		inline
558	IND	**MersenneTwisterRandomVariateGenerator::MersenneTwisterRandomVariateGenerator()
559	IND	**{
560	IND	****SetSeed ( 121212 );
561	IND	**}
562
563
564		/* Change log from MTRand.h */
565		// Change log:
566		//
567		// v0.1 - First release on 15 May 2000
568		//      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
569		//      - Translated from C to C++
570		//      - Made completely ANSI compliant
571		//      - Designed convenient interface for initialization, seeding, and
572		//        obtaining numbers in default or user-defined ranges
573		//      - Added automatic seeding from /dev/urandom or time() and clock()
574		//      - Provided functions for saving and loading generator state
575		//
576		// v0.2 - Fixed bug which reloaded generator one step too late
577		//
578		// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
579		//
580		// v0.4 - Removed trailing newline in saved generator format to be consistent
581		//        with output format of built-in types
582		//
583		// v0.5 - Improved portability by replacing static const int's with enum's and
584		//        clarifying return values in seed(); suggested by Eric Heimburg
585		//      - Removed MAXINT constant; use 0xffffffffUL instead
586		//
587		// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
588		//      - Changed integer [0,n] generator to give better uniformity
589		//
590		// v0.7 - Fixed operator precedence ambiguity in reload()
591		//      - Added access for real numbers in (0,1) and (0,n)
592		//
593		// v0.8 - Included time.h header to properly support time_t and clock_t
594		//
595		// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
596		//      - Allowed for seeding with arrays of any length
597		//      - Added access for real numbers in [0,1) with 53-bit resolution
598		//      - Added access for real numbers from normal (Gaussian) distributions
599		//      - Increased overall speed by optimizing twist()
600		//      - Doubled speed of integer [0,n] generation
601		//      - Fixed out-of-range number generation on 64-bit machines
602		//      - Improved portability by substituting literal constants for long enum's
603		//      - Changed license from GNU LGPL to BSD
604
605		} // end of namespace Statistics
606		} // end of namespace itk
607
608		#endif
609
610	EOF
611	EOF,EML

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:21PM

© Kitware Inc.