.. _program_listing_file_src_utils_RandomDist.hpp:

Program Listing for File RandomDist.hpp
=======================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_utils_RandomDist.hpp>` (``src/utils/RandomDist.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef RANDOMDIST_HPP
   #define RANDOMDIST_HPP
   
   #include "RngThread.hpp"
   #include <algorithm>
   #include <cmath>
   #include <cstddef>
   #include <cstdint>
   #include <limits>
   #include <numeric>
   #include <span>
   #include <stdexcept>
   #include <vector>
   
   #if defined(__GNUC__) && !defined(__clang__)
   #   pragma GCC push_options
   #   pragma GCC optimize("fp-contract=off")
   #endif
   
   namespace utils {
   
       inline auto canonical(RngType& rng) -> double
       {
           static_assert(RngType::min() == 0 &&
                         RngType::max() == std::numeric_limits<std::uint64_t>::max(),
                         "canonical() assumes an engine with a full 64-bit range");
           constexpr double twoTo64 = 0x1p64;
           const double u = static_cast<double>(rng()) / twoTo64;
           if (u >= 1.0) [[unlikely]] { return std::nextafter(1.0, 0.0); }
           return u;
       }
   
       class NormalDist {
   
       public:
   
           explicit NormalDist(const double mean = 0.0, const double stddev = 1.0)
               : mean_(mean), stddev_(stddev) {}
   
           auto operator()(RngType& rng) -> double
           {
   #ifdef __clang__
   #   pragma clang fp contract(off)
   #endif
               double ret = NAN;
               if (savedAvailable_)
               {
                   savedAvailable_ = false;
                   ret = saved_;
               }
               else
               {
                   double x = NAN;
                   double y = NAN;
                   double r2 = NAN;
                   while (true)
                   {
                       x = (2.0 * canonical(rng)) - 1.0;
                       y = (2.0 * canonical(rng)) - 1.0;
                       r2 = (x * x) + (y * y);
                       if (r2 <= 1.0 && r2 != 0.0) { break; }
                   }
                   const double mult = std::sqrt(-2.0 * std::log(r2) / r2);
                   saved_ = x * mult;
                   savedAvailable_ = true;
                   ret = y * mult;
               }
               return (ret * stddev_) + mean_;
           }
   
       private:
   
           double mean_;                  
           double stddev_;                
           double saved_ = 0.0;           
           bool savedAvailable_ = false;  
       };
   
       class LognormalDist {
   
       public:
   
           explicit LognormalDist(const double m = 0.0, const double s = 1.0)
               : m_(m), s_(s) {}
   
           auto operator()(RngType& rng) -> double
           {
   #ifdef __clang__
   #   pragma clang fp contract(off)
   #endif
               return std::exp((s_ * nd_(rng)) + m_);
           }
   
       private:
   
           double m_;       
           double s_;       
           NormalDist nd_;  
       };
   
       class DiscreteDist {
   
       public:
   
           explicit DiscreteDist(const std::span<const double> weights)
           {
               if (weights.size() < 2) { return; }
               const double sum = std::accumulate(weights.begin(), weights.end(), 0.0);
               if (!std::isfinite(sum) || sum <= 0.0)
               {
                   throw std::invalid_argument(
                       "DiscreteDist: sum of weights must be positive and finite");
               }
               cdf_.resize(weights.size());
               double cum = 0.0;
               for (std::size_t i = 0; i < weights.size(); ++i)
               {
                   cum += weights[i] / sum; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
                   cdf_[i] = cum;           // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
               }
               cdf_.back() = 1.0;
           }
   
           auto operator()(RngType& rng) const -> std::size_t
           {
               if (cdf_.empty()) { return 0; }
               const double p = canonical(rng);
               return static_cast<std::size_t>(
                   std::lower_bound(cdf_.begin(), cdf_.end(), p) - cdf_.begin());
           }
   
       private:
   
           std::vector<double> cdf_;  
       };
   
       class PoissonDist {
   
       public:
   
           explicit PoissonDist(const double mean) : mean_(mean)
           {
   #ifdef __clang__
   #   pragma clang fp contract(off)
   #endif
               if (!std::isfinite(mean) || mean < 0.0)
               {
                   throw std::invalid_argument(
                       "PoissonDist: mean must be non-negative and finite");
               }
               if (mean_ >= largeMean)
               {
                   const double m = std::floor(mean_);
                   lmThr_ = std::log(mean_);
                   lfm_ = std::lgamma(m + 1.0); // NOLINT(concurrency-mt-unsafe) -- only races on the global signgam, which is never read; libstdc++'s std::poisson_distribution calls lgamma identically
                   sm_ = std::sqrt(m);
                   const double dx = std::sqrt(2.0 * m * std::log(32.0 * m / piOver4));
                   d_ = std::round(std::max(6.0, std::min(m, dx)));
                   const double cx = (2.0 * m) + d_;
                   scx_ = std::sqrt(cx / 2.0);
                   oneOverCx_ = 1.0 / cx;
                   c2b_ = std::sqrt(piOver4 * cx) * std::exp(oneOverCx_);
                   cb_ = 2.0 * cx * std::exp(-d_ * oneOverCx_ * (1.0 + (d_ / 2.0))) / d_;
               }
               else
               {
                   lmThr_ = std::exp(-mean_);
               }
           }
   
           auto operator()(RngType& rng) -> unsigned int
           {
   #ifdef __clang__
   #   pragma clang fp contract(off)
   #endif
               if (mean_ < largeMean)
               {
                   unsigned int x = 0;
                   double prod = 1.0;
                   while (true)
                   {
                       prod *= canonical(rng);
                       x += 1;
                       if (prod <= lmThr_) { break; }
                   }
                   return x - 1;
               }
   
               // Large mean: rejection method
               constexpr double naf = (1.0 - std::numeric_limits<double>::epsilon()) / 2.0;
               constexpr double thr =
                   static_cast<double>(std::numeric_limits<unsigned int>::max()) + naf;
               constexpr double sqrtPiOver2 = 1.2533141373155002512;
               constexpr double oneOver78 = 1.0 / 78.0;
               constexpr double expOneOver78 = 1.0129030479320018583;
               const double m = std::floor(mean_);
               const double c1 = sm_ * sqrtPiOver2;
               const double c2 = c2b_ + c1;
               const double c3 = c2 + 1.0;
               const double c4 = c3 + 1.0;
               const double c5 = c4 + expOneOver78;
               const double c = cb_ + c5;
               const double twoCx = 2.0 * ((2.0 * m) + d_);
   
               while (true)
               {
                   const double u = c * canonical(rng);
                   const double e = -std::log(1.0 - canonical(rng));
                   double x = 0.0;
                   double w = 0.0;
                   if (u <= c1)
                   {
                       const double n = nd_(rng);
                       const double y = (-std::abs(n) * sm_) - 1.0;
                       x = std::floor(y);
                       w = -n * n / 2.0;
                       if (x < -m) { continue; }
                   }
                   else if (u <= c2)
                   {
                       const double n = nd_(rng);
                       const double y = 1.0 + (std::abs(n) * scx_);
                       x = std::ceil(y);
                       w = y * (2.0 - y) * oneOverCx_;
                       if (x > d_) { continue; }
                   }
                   else if (u <= c3)
                   {
                       x = -1.0;
                   }
                   else if (u <= c4)
                   {
                       x = 0.0;
                   }
                   else if (u <= c5)
                   {
                       x = 1.0;
                       w = oneOver78;
                   }
                   else
                   {
                       const double v = -std::log(1.0 - canonical(rng));
                       const double y = d_ + (v * twoCx / d_);
                       x = std::ceil(y);
                       w = -d_ * oneOverCx_ * (1.0 + (y / 2.0));
                   }
                   const bool reject =
                       (w - e - (x * lmThr_) > lfm_ - std::lgamma(x + m + 1.0)) || // NOLINT(concurrency-mt-unsafe) -- see constructor
                       (x + m >= thr);
                   if (!reject)
                   {
                       return static_cast<unsigned int>(x + m + naf);
                   }
               }
           }
   
       private:
   
           static constexpr double largeMean = 12.0;  
           static constexpr double piOver4 = 0.78539816339744830962;  
           double mean_;             
           double lmThr_ = 0.0;      
           double lfm_ = 0.0;        
           double sm_ = 0.0;         
           double d_ = 0.0;          
           double scx_ = 0.0;        
           double oneOverCx_ = 0.0;  
           double c2b_ = 0.0;        
           double cb_ = 0.0;         
           NormalDist nd_;           
       };
   
   } // namespace utils
   
   #if defined(__GNUC__) && !defined(__clang__)
   #   pragma GCC pop_options
   #endif
   
   #endif // RANDOMDIST_HPP