Program Listing for File RandomDist.hpp
↰ Return to documentation for file (src/utils/RandomDist.hpp)
#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