itkMersenneTwisterRandomVariateGenerator.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkMersenneTwisterRandomVariateGenerator_h
#define itkMersenneTwisterRandomVariateGenerator_h
 
#include "itkMacro.h"
#include "itkObjectFactory.h"
#include "itkRandomVariateGeneratorBase.h"
#include "itkIntTypes.h"
#include "itkMath.h"
#include "itkSingletonMacro.h"
 
#include <atomic>
#include <mutex>
#include <climits>
#include <ctime>
#include <limits>
 
namespace itk
{
namespace Statistics
{
 
struct MersenneTwisterGlobals;
 
class ITKCommon_EXPORT MersenneTwisterRandomVariateGenerator : public RandomVariateGeneratorBase
{
public:
  using Self = MersenneTwisterRandomVariateGenerator;
  using Superclass = RandomVariateGeneratorBase;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
 
  using IntegerType = uint32_t;

  itkOverrideGetNameOfClassMacro(MersenneTwisterRandomVariateGenerator);

  static Pointer
  New();

  static Pointer
  GetInstance();

  static void
  ResetNextSeed();

  static constexpr IntegerType StateVectorLength = 624;

  void
  Initialize(const IntegerType seed);
 
  /* Initialize with clock time */
  void
  Initialize()
  {
    SetSeed();
  }

  double
  GetVariateWithClosedRange()
  {
    return static_cast<double>(GetIntegerVariate()) * (1.0 / double{ std::numeric_limits<uint32_t>::max() });
  }

  double
  GetVariateWithClosedRange(const double n)
  {
    return GetVariateWithClosedRange() * n;
  }

  double
  GetVariateWithOpenUpperRange()
  {
    return static_cast<double>(GetIntegerVariate()) / double{ 1ULL << 32ULL };
  }

  double
  GetVariateWithOpenUpperRange(const double n)
  {
    return GetVariateWithOpenUpperRange() * n;
  }

  double
  GetVariateWithOpenRange()
  {
    return (static_cast<double>(GetIntegerVariate()) + 0.5) / double{ 1ULL << 32ULL };
  }

  double
  GetVariateWithOpenRange(const double n)
  {
    return GetVariateWithOpenRange() * n;
  }

  IntegerType
  GetIntegerVariate();

  IntegerType
  GetIntegerVariate(const IntegerType n);

  double
  Get53BitVariate();
 
  /* Access to a normal random number distribution
   * TODO: Compare with vnl_sample_normal */
  double
  GetNormalVariate(const double mean = 0.0, const double variance = 1.0);
 
  /* Access to a uniform random number distribution in the range [a, b)
   * TODO: Compare with vnl_sample_uniform */
  double
  GetUniformVariate(const double a, const double b);

  double
  GetVariate() override
  {
    return GetVariateWithClosedRange();
  }

  double
  operator()()
  {
    return GetVariate();
  }

  void
  SetSeed(const IntegerType oneSeed)
  {
    // Seed the generator with a simple IntegerType
    Initialize(oneSeed);
  }
 
  void
  SetSeed()
  {
    SetSeed(Self::CreateRandomSeed());
  }

  IntegerType
  GetSeed() const
  {
    return m_Seed;
  }

  static IntegerType
  GetNextSeed();
 
  /*
  // Saving and loading generator state
  void save( IntegerType* saveArray ) const;  // to array of size SAVE
  void load( IntegerType *const loadArray );  // from such array
  */
 
protected:
  MersenneTwisterRandomVariateGenerator();
  ~MersenneTwisterRandomVariateGenerator() override;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;

  void
  reload();
 
  static IntegerType
  hash(time_t t, clock_t c);
 
private:
  itkGetGlobalDeclarationMacro(MersenneTwisterGlobals, PimplGlobals);

  static Pointer
  CreateInstance();

  static IntegerType
  CreateRandomSeed()
  {
    return hash(time(nullptr), clock());
  }

  void
  InitializeWithoutMutexLocking(const IntegerType seed);
 
  // Internal state
  IntegerType m_State[StateVectorLength];
 
  // Next value to get from state
  IntegerType * m_PNext{};
 
  // Number of values left before reload is needed
  int m_Left{};
 
  // Seed value
  std::atomic<IntegerType> m_Seed{};
 
  // Local lock to enable concurrent access to singleton
  std::mutex m_InstanceMutex{};
 
  // Static/Global Variable need to be thread-safely accessed
 
  static MersenneTwisterGlobals * m_PimplGlobals;
 
}; // end of class
 
// Declare inlined functions.... (must be declared in the header)

inline double
MersenneTwisterRandomVariateGenerator::Get53BitVariate()
{
  const IntegerType a = GetIntegerVariate() >> 5;
  const IntegerType b = GetIntegerVariate() >> 6;
 
  return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku
                                                            // Wada
}

// TODO: Compare with vnl_sample_normal
inline double
MersenneTwisterRandomVariateGenerator::GetNormalVariate(const double mean, const double variance)
{
  // Return a real number from a normal (Gaussian) distribution with given
  // mean and variance by Box-Muller method
  const double r = std::sqrt(-2.0 * std::log(1.0 - GetVariateWithOpenRange()) * variance);
  const double phi = 2.0 * itk::Math::pi * GetVariateWithOpenUpperRange();
 
  return mean + r * std::cos(phi);
}

// TODO: Compare with vnl_sample_uniform
inline double
MersenneTwisterRandomVariateGenerator::GetUniformVariate(const double a, const double b)
{
  const double u = GetVariateWithOpenUpperRange();
 
  return ((1.0 - u) * a + u * b);
}
 
/* Change log from MTRand.h */
// Change log:
//
// v0.1 - First release on 15 May 2000
//      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
//      - Translated from C to C++
//      - Made completely ANSI compliant
//      - Designed convenient interface for initialization, seeding, and
//        obtaining numbers in default or user-defined ranges
//      - Added automatic seeding from /dev/urandom or time() and clock()
//      - Provided functions for saving and loading generator state
//
// v0.2 - Fixed bug which reloaded generator one step too late
//
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
//
// v0.4 - Removed trailing newline in saved generator format to be consistent
//        with output format of built-in types
//
// v0.5 - Improved portability by replacing static const int's with enum's and
//        clarifying return values in seed(); suggested by Eric Heimburg
//      - Removed MAXINT constant; use 0xffffffffUL instead
//
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
//      - Changed integer [0,n] generator to give better uniformity
//
// v0.7 - Fixed operator precedence ambiguity in reload()
//      - Added access for real numbers in (0,1) and (0,n)
//
// v0.8 - Included time.h header to properly support time_t and clock_t
//
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
//      - Allowed for seeding with arrays of any length
//      - Added access for real numbers in [0,1) with 53-bit resolution
//      - Added access for real numbers from normal (Gaussian) distributions
//      - Increased overall speed by optimizing twist()
//      - Doubled speed of integer [0,n] generation
//      - Fixed out-of-range number generation on 64-bit machines
//      - Improved portability by substituting literal constants for long enum's
//      - Changed license from GNU LGPL to BSD
} // end namespace Statistics
} // end namespace itk
 
#endif