ITK  6.0.0
Insight Toolkit
itkMersenneTwisterRandomVariateGenerator.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright NumFOCUS
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 * https://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
33namespace itk
34{
35namespace Statistics
36{
125struct MersenneTwisterGlobals;
126
128{
129public:
135
136 using IntegerType = uint32_t;
137
139 itkOverrideGetNameOfClassMacro(MersenneTwisterRandomVariateGenerator);
140
149 static Pointer
151
160 static Pointer
162
165 static void
167
169 static constexpr IntegerType StateVectorLength = 624;
170
172 void
173 Initialize(const IntegerType seed);
174
175 /* Initialize with clock time */
176 void
177 Initialize();
178
180 double
181 GetVariateWithClosedRange();
182
184 double
185 GetVariateWithClosedRange(const double n);
186
188 double
189 GetVariateWithOpenUpperRange();
190
192 double
193 GetVariateWithOpenUpperRange(const double n);
194
196 double
197 GetVariateWithOpenRange();
198
200 double
201 GetVariateWithOpenRange(const double n);
202
205 GetIntegerVariate();
206
209 GetIntegerVariate(const IntegerType & n);
210
213 double
214 Get53BitVariate();
215
216 /* Access to a normal random number distribution
217 * TODO: Compare with vnl_sample_normal */
218 double
219 GetNormalVariate(const double mean = 0.0, const double variance = 1.0);
220
221 /* Access to a uniform random number distribution in the range [a, b)
222 * TODO: Compare with vnl_sample_uniform */
223 double
224 GetUniformVariate(const double a, const double b);
225
231 double
232 GetVariate() override;
233
235 double
236 operator()();
237
242 inline void
243 SetSeed(const IntegerType oneSeed);
244 inline void
245 SetSeed();
253 GetSeed() const;
254
260 static IntegerType
262
263 /*
264 // Saving and loading generator state
265 void save( IntegerType* saveArray ) const; // to array of size SAVE
266 void load( IntegerType *const loadArray ); // from such array
267 */
268
269protected:
272 void
273 PrintSelf(std::ostream & os, Indent indent) const override;
274
276 static constexpr unsigned int M = 397;
277
279 void
280 reload();
281
283 hiBit(const IntegerType & u) const
284 {
285 return u & 0x80000000;
286 }
287 IntegerType
288 loBit(const IntegerType & u) const
289 {
290 return u & 0x00000001;
291 }
292 IntegerType
293 loBits(const IntegerType & u) const
294 {
295 return u & 0x7fffffff;
296 }
297 IntegerType
298 mixBits(const IntegerType & u, const IntegerType & v) const
299 {
300 return hiBit(u) | loBits(v);
301 }
302
303 IntegerType
304 twist(const IntegerType & m, const IntegerType & s0, const IntegerType & s1) const
305 {
306 return m ^ (mixBits(s0, s1) >> 1) ^ (-static_cast<int32_t>(loBit(s1)) & 0x9908b0df);
307 }
308
309 static IntegerType
310 hash(time_t t, clock_t c);
311
312 // Internal state
313 IntegerType state[StateVectorLength];
314
315 // Next value to get from state
316 IntegerType * m_PNext{};
317
318 // Number of values left before reload is needed
319 int m_Left{};
320
321 // Seed value
322 std::atomic<IntegerType> m_Seed{};
323
324private:
326 itkGetGlobalDeclarationMacro(MersenneTwisterGlobals, PimplGlobals);
327
329 static Pointer
331
332 // Local lock to enable concurrent access to singleton
333 std::mutex m_InstanceMutex{};
334
335 // Static/Global Variable need to be thread-safely accessed
336
337 static MersenneTwisterGlobals * m_PimplGlobals;
338
339}; // end of class
340
341// Declare inlined functions.... (must be declared in the header)
342
343inline void
345{
346 const std::lock_guard<std::mutex> lockGuard(m_InstanceMutex);
347 this->m_Seed = seed;
348 // Initialize generator state with seed
349 // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
350 // In previous versions, most significant bits (MSBs) of the seed affect
351 // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
352 IntegerType * s = state;
353 IntegerType * r = state;
354 IntegerType i = 1;
355
356 *s++ = seed & 0xffffffffUL;
358 {
359 *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
360 ++r;
361 }
362 reload();
363}
364
365inline void
367{
368 // Generate N new values in state
369 // Made clearer and faster by Matthew Bellew
370 // matthew dot bellew at home dot com
371
372 // get rid of VS warning
373 constexpr auto index = int{ M } - int{ MersenneTwisterRandomVariateGenerator::StateVectorLength };
374
375 IntegerType * p = state;
376 int i;
377
379 {
380 *p = twist(p[M], p[0], p[1]);
381 }
382 for (i = M; --i; ++p)
383 {
384 *p = twist(p[index], p[0], p[1]);
385 }
386 *p = twist(p[index], p[0], state[0]);
387
389 m_PNext = state;
390}
391
392inline void
394{
395 SetSeed();
396}
397
398inline void
400{
401 // Seed the generator with a simple IntegerType
402 Initialize(oneSeed);
403}
404
405inline void
407{
408 // use time() and clock() to generate a unlikely-to-repeat seed.
409 SetSeed(hash(time(nullptr), clock()));
410}
411
412
415{
416 return this->m_Seed;
417}
418
422{
423 if (m_Left == 0)
424 {
425 reload();
426 }
427 --m_Left;
430 IntegerType s1 = *m_PNext++;
431 s1 ^= (s1 >> 11);
432 s1 ^= (s1 << 7) & 0x9d2c5680;
433 s1 ^= (s1 << 15) & 0xefc60000;
434 return (s1 ^ (s1 >> 18));
435}
436
437inline double
439{
440 return static_cast<double>(GetIntegerVariate()) * (1.0 / 4294967295.0);
441}
442
444inline double
446{
447 return GetVariateWithClosedRange() * n;
448}
449
451inline double
453{
454 return static_cast<double>(GetIntegerVariate()) * (1.0 / 4294967296.0);
455}
456
458inline double
460{
461 return GetVariateWithOpenUpperRange() * n;
462}
463
465inline double
467{
468 return (static_cast<double>(GetIntegerVariate()) + 0.5) * (1.0 / 4294967296.0);
469}
470
472inline double
474{
475 return GetVariateWithOpenRange() * n;
476}
477
480{
481 // Find which bits are used in n
482 IntegerType used = n;
483
484 used |= used >> 1;
485 used |= used >> 2;
486 used |= used >> 4;
487 used |= used >> 8;
488 used |= used >> 16;
489
490 // Draw numbers until one is found in [0,n]
491 IntegerType i;
492 do
493 {
494 i = GetIntegerVariate() & used; // toss unused bits to shorten search
495 } while (i > n);
496
497 return i;
498}
499
502inline double
504{
505 IntegerType a = GetIntegerVariate() >> 5, b = GetIntegerVariate() >> 6;
506
507 return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku
508 // Wada
509}
510
512// TODO: Compare with vnl_sample_normal
513inline double
514MersenneTwisterRandomVariateGenerator::GetNormalVariate(const double mean, const double variance)
515{
516 // Return a real number from a normal (Gaussian) distribution with given
517 // mean and variance by Box-Muller method
518 double r = std::sqrt(-2.0 * std::log(1.0 - GetVariateWithOpenRange()) * variance);
519 double phi = 2.0 * itk::Math::pi * GetVariateWithOpenUpperRange();
522 return mean + r * std::cos(phi);
523}
524
526// TODO: Compare with vnl_sample_uniform
527inline double
529{
530 double u = GetVariateWithOpenUpperRange();
531
532 return ((1.0 - u) * a + u * b);
533}
534
535inline double
537{
539}
540
541inline double
543{
544 return GetVariate();
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 namespace Statistics
588} // end namespace itk
589
590#endif
Control indentation during Print() invocation.
Definition: itkIndent.h:50
Light weight base class for most itk classes.
IntegerType mixBits(const IntegerType &u, const IntegerType &v) const
static Pointer New()
Method for creation through the object factory.
void PrintSelf(std::ostream &os, Indent indent) const override
IntegerType twist(const IntegerType &m, const IntegerType &s0, const IntegerType &s1) const
static IntegerType hash(time_t t, clock_t c)
itkGetGlobalDeclarationMacro(MersenneTwisterGlobals, PimplGlobals)
double GetNormalVariate(const double mean=0.0, const double variance=1.0)
Defines common interfaces for random variate generators.
static constexpr double pi
Definition: itkMath.h:66
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....