Random.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/Random.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_Random_h
#define __PCL_Random_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/AutoPointer.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
extern uint64 RandomSeed64();
 
inline uint32 RandomSeed32()
{
   union { uint64 u64; uint32 u32[ 2 ]; } seed;
   seed.u64 = RandomSeed64();
   return seed.u32[0] ^ seed.u32[1];
}
 
// ----------------------------------------------------------------------------
 
class FastMersenneTwister;
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS RandomNumberGenerator
{
public:
 
   RandomNumberGenerator( double ymax = 1.0, uint32 seed = 0 );
 
   virtual ~RandomNumberGenerator();
 
   double operator ()()
   {
      return m_rmax*Rand32();
   }
 
   uint32 Rand32();
 
   double Rand1()
   {
      return double( Rand32() )/uint32_max;
   }
 
   double Uniform()
   {
      return operator()();
   }
 
   double Normal( double mean = 0, double sigma = 1 );
 
   double Gaussian( double mean = 0, double sigma = 1 )
   {
      return Normal( mean, sigma );
   }
 
   int Poisson( double lambda );
 
   double UpperBound() const
   {
      return m_ymax;
   }
 
   void SetUpperBound( double ymax )
   {
      PCL_PRECONDITION( ymax > 0 )
      PCL_PRECONDITION( 1 + ymax != 1 )
      m_rmax = (m_ymax = ymax)/double( uint32_max );
      m_normal = false;
   }
 
private:
 
   AutoPointer<FastMersenneTwister> m_generator;
   double                           m_ymax;
   double                           m_rmax;
   bool                             m_normal;
   double                           m_vs;     // second result from Box–Muller transform
   DVector                          m_lambda; // precalculated for current Poisson lambda
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS XorShift1024
{
public:
 
   XorShift1024( uint64 seed = 0 ) noexcept( false )
   {
      Initialize( seed );
   }
 
   double operator ()() noexcept
   {
      return 5.4210108624275221703311e-20 * UI64(); // 1.0/(2^64 -1)
   }
 
   uint64 UI64() noexcept
   {
      uint64 s0 = m_s[m_p];
      uint64 s1 = m_s[m_p = (m_p + 1) & 15];
      s1 ^= s1 << 31; // a
      s1 ^= s1 >> 11; // b
      s0 ^= s0 >> 30; // c
      return (m_s[m_p] = s0 ^ s1) * 1181783497276652981ull;
   }
 
   uint32 UI32() noexcept
   {
      return uint32( UI64() );
   }
 
   uint64 UI64N( uint64 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UIN( uint32 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UI32N( uint32 n ) noexcept
   {
      return UIN( n );
   }
 
   void Initialize( uint64 x )
   {
      if ( x == 0 )
         x = RandomSeed64();
      // Use a xorshift64* generator to initialize the state space.
      for ( int i = 0; i < 16; ++i )
      {
         x ^= x >> 12; // a
         x ^= x << 25; // b
         x ^= x >> 27; // c
         m_s[i] = x * 2685821657736338717ull;
      }
      m_p = 0;
   }
 
private:
 
   uint64 m_s[ 16 ]; // state space
   int m_p;          // current index
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS XoShiRoBase
{
public:
 
   XoShiRoBase() = default;
 
protected:
 
   static uint64 RotL( const uint64 x, int k ) noexcept
   {
      return (x << k) | (x >> (64 - k));
   }
 
   static uint64 SplitMix64( uint64& x ) noexcept
   {
      uint64 z = (x += 0x9e3779b97f4a7c15);
      z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
      z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
      return z ^ (z >> 31);
   }
 
   static double UI64ToDouble( uint64 x ) noexcept
   {
      return (x >> 11) * 0x1.0p-53;
   }
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS XoShiRo256ss : public XoShiRoBase
{
public:
 
   XoShiRo256ss( uint64 seed = 0 ) noexcept( false )
   {
      Initialize( seed );
   }
 
   double operator ()() noexcept
   {
      return UI64ToDouble( UI64() );
   }
 
   uint64 UI64() noexcept
   {
      const uint64 result = RotL( m_s[1]*5, 7 ) * 9;
      const uint64 t = m_s[1] << 17;
      m_s[2] ^= m_s[0];
      m_s[3] ^= m_s[1];
      m_s[1] ^= m_s[2];
      m_s[0] ^= m_s[3];
      m_s[2] ^= t;
      m_s[3] = RotL( m_s[3], 45 );
      return result;
   }
 
   uint32 UI32() noexcept
   {
      return uint32( UI64() );
   }
 
   uint64 UI64N( uint64 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UIN( uint32 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UI32N( uint32 n ) noexcept
   {
      return UIN( n );
   }
 
   void Initialize( uint64 x )
   {
      if ( x == 0 )
         x = RandomSeed64();
      // Use a SplitMix64 generator to initialize the state space.
      for ( int i = 0; i < 4; ++i )
         m_s[i] = SplitMix64( x );
   }
 
private:
 
   uint64 m_s[ 4 ];
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS XoRoShiRo1024ss : public XoShiRoBase
{
public:
 
   XoRoShiRo1024ss( uint64 seed = 0 ) noexcept( false )
   {
      Initialize( seed );
   }
 
   double operator ()() noexcept
   {
      return UI64ToDouble( UI64() );
   }
 
   uint64 UI64() noexcept
   {
      const int q = m_p;
      const uint64 s0 = m_s[m_p = (m_p + 1) & 15];
      uint64 s15 = m_s[q];
      const uint64 result = RotL( s0*5, 7 ) * 9;
      s15 ^= s0;
      m_s[q] = RotL( s0, 25 ) ^ s15 ^ (s15 << 27);
      m_s[m_p] = RotL( s15, 36 );
      return result;
   }
 
   uint32 UI32() noexcept
   {
      return uint32( UI64() );
   }
 
   uint64 UI64N( uint64 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UIN( uint32 n ) noexcept
   {
      return UI64() % n;
   }
 
   uint32 UI32N( uint32 n ) noexcept
   {
      return UIN( n );
   }
 
   void Initialize( uint64 x )
   {
      if ( x == 0 )
         x = RandomSeed64();
      // Use a SplitMix64 generator to initialize the state space.
      for ( int i = 0; i < 16; ++i )
         m_s[i] = SplitMix64( x );
      m_p = 0;
   }
 
private:
 
   uint64 m_s[ 16 ];
   int    m_p;
};
 
// ----------------------------------------------------------------------------
 
template <class RNG>
class PCL_CLASS NormalRandomDeviates
{
public:
 
   NormalRandomDeviates( RNG& R ) noexcept( false )
      : m_R( R )
   {
   }
 
   double operator ()() noexcept
   {
      /*
       * Marsaglia polar method.
       */
      double x;
      if ( m_first )
      {
         do
         {
            double u1 = m_R();
            double u2 = m_R();
            m_v1 = 2*u1 - 1;
            m_v2 = 2*u2 - 1;
            m_s = m_v1*m_v1 + m_v2*m_v2;
         }
         while ( m_s >= 1 || m_s <= std::numeric_limits<double>::epsilon() );
 
         x = m_v1 * Sqrt( -2*Ln( m_s )/m_s );
      }
      else
         x = m_v2 * Sqrt( -2*Ln( m_s )/m_s );
 
      m_first = !m_first;
      return x;
   }
 
private:
 
   RNG&   m_R;
   double m_v1 = 0;
   double m_v2 = 0;
   double m_s = 0;
   bool   m_first = true;
};
 
// ----------------------------------------------------------------------------
 
template <class RNG>
class PCL_CLASS PoissonRandomDeviates
{
public:
 
   PoissonRandomDeviates( RNG& R ) noexcept( false )
      : m_R( R )
   {
   }
 
   int operator ()( double value ) noexcept
   {
      if ( value < 30 )
      {
         /*
          * Implementation of the algorithm by Donald E. Knuth, 1969.
          *
          * This algorithm is slow (unusable) for large values.
          */
         double p = 1, L = Exp( -value );
         int k = 0;
         do
         {
            ++k;
            p *= m_R();
         }
         while ( p > L );
         return k-1;
      }
 
      /*
       * Code adapted from 'Random number generation in C++', by John D. Cook:
       *
       * https://www.johndcook.com/blog/cpp_random_number_generation/
       *
       * The algorithm is from "The Computer Generation of Poisson Random
       * Variables" by A. C. Atkinson, Journal of the Royal Statistical
       * Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
       *
       * This algorithm is slow (unusable) for small values.
       */
      double c = 0.767 - 3.36/value;
      double beta = Const<double>::pi()/Sqrt( 3*value );
      double alpha = beta*value;
      double k = Ln( c ) - value - Ln( beta );
      for ( ;; )
      {
         double u = m_R();
         double x = (alpha - Ln( (1 - u)/u ))/beta;
         int n = int( Floor( x + 0.5 ) );
         if ( n < 0 )
            continue;
         double v = m_R();
         double y = alpha - beta*x;
         double temp = 1 + Exp( y );
         double lhs = y + Ln( v/temp/temp );
         double rhs = k + n*Ln( value ) - LnFactorial( n );
         if ( lhs <= rhs )
            return n;
      }
   }
 
private:
 
   RNG& m_R;
};
 
// ----------------------------------------------------------------------------
 
template <class RNG>
class PCL_CLASS GammaRandomDeviates
{
public:
 
   GammaRandomDeviates( RNG& R, double shape = 1, double scale = 1 ) noexcept( false )
      : m_R( R )
      , m_shape( shape )
      , m_scale( scale )
      , m_normal( R )
   {
      if ( m_shape <= 0 )
         throw Error( "GammaRandomDeviates(): The function shape parameter must be > 0." );
      if ( m_scale <= 0 )
         throw Error( "GammaRandomDeviates(): The scale parameter must be > 0." );
 
      m_d = ((m_shape >= 1) ? m_shape : m_shape + 1) - 1.0/3.0;
      m_c = 1/Sqrt( 9*m_d );
   }
 
   double operator ()() noexcept
   {
      /*
       * Code adapted from 'Random number generation in C++', by John D. Cook:
       *    https://www.johndcook.com/blog/cpp_random_number_generation/
       *
       * Implementation based on "A Simple Method for Generating Gamma
       * Variables" by George Marsaglia and Wai Wan Tsang. ACM Transactions on
       * Mathematical Software Vol 26, No 3, September 2000, pages 363-372.
       */
      for ( ;; )
      {
         double x, v;
         do
         {
            x = m_normal();
            v = 1 + m_c*x;
         }
         while ( v <= 0 );
         v = v*v*v;
         double u = m_R();
         double xsquared = x*x;
         if ( u < 1 - 0.0331*xsquared*xsquared || Ln( u ) < 0.5*xsquared + m_d*(1 - v + Ln( v )) )
         {
            double g = m_scale*m_d*v;
            if ( m_shape < 1 )
               g *= Pow( m_R(), 1/m_shape );
            return g;
         }
      }
   }
 
private:
 
   RNG&                      m_R;
   double                    m_shape;
   double                    m_scale;
   double                    m_d;
   double                    m_c;
   NormalRandomDeviates<RNG> m_normal;
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_Random_h
 
// ----------------------------------------------------------------------------
// EOF pcl/Random.h - Released 2024-06-18T15:48:54Z