Utilities/detail/FpeSafeGammaDistribution.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cmath>
0012 #include <limits>
0013 #include <random>
0014
0015 namespace ActsFatras::detail {
0016
0017 /// Draw random numbers from a gamma distribution.
0018 /// Does not generate FPE underflow
0019 ///
0020 /// Implements the same interface as the standard library distributions.
0021 class FpeSafeGammaDistribution {
0022  public:
0023   /// Parameter struct that contains all distribution parameters.
0024   struct param_type {
0025     /// Parameters must link back to the host distribution.
0026     using distribution_type = FpeSafeGammaDistribution;
0027
0028     /// Location parameter.
0029     ///
0030     /// Shape parameter
0031     double alpha = 0.5;
0032     /// Scale parameter.
0033     double theta = 1.0;
0034
0035     /// Construct from parameters.
0036     param_type(double alpha_, double theta_) : alpha(alpha_), theta(theta_) {}
0037     // Explicitly defaulted construction and assignment
0038     param_type() = default;
0039     param_type(const param_type &) = default;
0040     param_type(param_type &&) = default;
0041     param_type &operator=(const param_type &) = default;
0042     param_type &operator=(param_type &&) = default;
0043
0044     /// Parameters should be EqualityComparable
0045     friend bool operator==(const param_type &lhs, const param_type &rhs) {
0046       return (lhs.alpha == rhs.alpha) && (lhs.theta == rhs.theta);
0047     }
0048   };
0049   /// The type of the generated values.
0050   using result_type = double;
0051
0052   /// Construct directly from the distribution parameters.
0053   FpeSafeGammaDistribution(double alpha, double theta) : m_cfg(alpha, theta) {}
0054   /// Construct from a parameter object.
0055   explicit FpeSafeGammaDistribution(const param_type &cfg) : m_cfg(cfg) {}
0056   // Explicitlely defaulted construction and assignment
0057   FpeSafeGammaDistribution() = default;
0058   FpeSafeGammaDistribution(const FpeSafeGammaDistribution &) = default;
0059   FpeSafeGammaDistribution(FpeSafeGammaDistribution &&) = default;
0060   FpeSafeGammaDistribution &operator=(const FpeSafeGammaDistribution &) =
0061       default;
0062   FpeSafeGammaDistribution &operator=(FpeSafeGammaDistribution &&) = default;
0063
0064   /// Reset any possible internal state. Noop, since there is no internal state.
0065   void reset() {}
0066   /// Return the currently configured distribution parameters.
0067   param_type param() const { return m_cfg; }
0068   /// Set the distribution parameters.
0069   void param(const param_type &cfg) { m_cfg = cfg; }
0070
0071   /// The minimum value the distribution generates.
0072   result_type min() const { return 0.; }
0073   /// The maximum value the distribution generates.
0074   result_type max() const { return std::numeric_limits<double>::infinity(); }
0075
0076   /// Generate a random number from the configured gamma distribution.
0077   template <typename Generator>
0078   result_type operator()(Generator &generator) {
0079     return (*this)(generator, m_cfg);
0080   }
0081   /// Generate a random number from the given gamma distribution.
0082   template <typename Generator>
0083   result_type operator()(Generator &generator, const param_type &params) {
0084     if (params.alpha >= 1.) {
0085       std::gamma_distribution<double> gDist(params.alpha, params.theta);
0086       return gDist(generator);
0087     } else if (params.alpha <= 0.) {
0088       return 0.;
0089     }
0090
0091     // This is from libstdc++; libcxx would give a different result
0092     std::gamma_distribution<double> gDistPlus1(params.alpha + 1., params.theta);
0093     // Get a random number using alpha+1
0094     const auto uPlus1 = gDistPlus1(generator);
0095     std::uniform_real_distribution<double> uDist(0., 1.);
0096     double u2 = uDist(generator);
0097     while (u2 == 0.) {
0098       u2 = uDist(generator);
0099     }
0100     // Check if this would cause underflow
0101     const double invAlpha = 1. / params.alpha;
0102     if (std::log(u2) * invAlpha + std::log(uPlus1) <
0103         std::log(std::numeric_limits<double>::min())) {
0104       // This would be underflow - get 0.0 right away
0105       return 0.;
0106     } else {
0107       return uPlus1 * std::pow(u2, invAlpha);
0108     }
0109   }
0110
0111   /// Provide standard comparison operators
0112   friend bool operator==(const FpeSafeGammaDistribution &lhs,
0113                          const FpeSafeGammaDistribution &rhs) {
0114     return lhs.m_cfg == rhs.m_cfg;
0115   }
0116
0117  private:
0118   param_type m_cfg;
0119 };
0120
0121 }  // namespace ActsFatras::detail