random/distribution/PoissonDistribution.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file corecel/random/distribution/PoissonDistribution.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "celeritas/Constants.hh"
0016
0017 #include "NormalDistribution.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Sample from a Poisson distribution.
0024  *
0025  * The Poisson distribution describes the probability of \f$ k \f$ events
0026  * occuring in a fixed interval given a mean rate of occurance \f$ \lambda \f$
0027  * and has the PMF:
0028  * \f[
0029    f(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!} \:.
0030    \f]
0031  * For small \f$ \lambda \f$, a direct method described in Knuth, Donald E.,
0032  * Seminumerical Algorithms, The Art of Computer Programming, Volume 2 can be
0033  * used to generate samples from the Poisson distribution. Uniformly
0034  * distributed random numbers are generated until the relation
0035  * \f[
0036    \prod_{k = 1}^n U_k \le e^{-\lambda}
0037    \f]
0038  * is satisfied; then, the random variable \f$ X = n - 1 \f$ will have a
0039  * Poisson distribution. On average this approach requires the generation of
0040  * \f$ \lambda + 1 \f$ uniform random samples, so a different method should be
0041  * used for large \f$ \lambda \f$.
0042  *
0043  * Geant4 uses Knuth's algorithm for \f$ \lambda \le 16 \f$ and a Gaussian
0044  * approximation for \f$ \lambda > 16 \f$ (see \c G4Poisson), which is faster
0045  * but less accurate than other methods. The same approach is used here.
0046  *
0047  * \todo Break this into two distributions: one actual poisson distribution,
0048  * one "integer normal" distribution, and a variant type that selects between
0049  * them. In most cases we care about, lambda is small.
0050  */
0051 template<class RealType = ::celeritas::real_type>
0052 class PoissonDistribution
0053 {
0054   public:
0055     //!@{
0056     //! \name Type aliases
0057     using real_type = RealType;
0058     using result_type = unsigned int;
0059     //!@}
0060
0061   public:
0062     // Construct with defaults
0063     explicit inline CELER_FUNCTION PoissonDistribution(real_type lambda = 1);
0064
0065     // Sample a random number according to the distribution
0066     template<class Generator>
0067     inline CELER_FUNCTION result_type operator()(Generator& rng);
0068
0069     //! Maximum value of lambda for using the direct method
0070     static CELER_CONSTEXPR_FUNCTION int lambda_threshold() { return 16; }
0071
0072   private:
0073     real_type const lambda_;
0074     NormalDistribution<real_type> sample_normal_;
0075 };
0076
0077 //---------------------------------------------------------------------------//
0078 // INLINE DEFINITIONS
0079 //---------------------------------------------------------------------------//
0080 /*!
0081  * Construct from the mean of the Poisson distribution.
0082  */
0083 template<class RealType>
0084 CELER_FUNCTION
0085 PoissonDistribution<RealType>::PoissonDistribution(real_type lambda)
0086     : lambda_(lambda), sample_normal_(lambda_, std::sqrt(lambda_))
0087 {
0088     CELER_EXPECT(lambda_ > 0);
0089 }
0090
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Sample a random number according to the distribution.
0094  */
0095 template<class RealType>
0096 template<class Generator>
0097 CELER_FUNCTION auto
0098 PoissonDistribution<RealType>::operator()(Generator& rng) -> result_type
0099 {
0100     if (lambda_ <= PoissonDistribution::lambda_threshold())
0101     {
0102         // Use direct method
0103         int k = 0;
0104         real_type p = std::exp(lambda_);
0105         do
0106         {
0107             ++k;
0108             p *= generate_canonical<real_type>(rng);
0109         } while (p > 1);
0110         return static_cast<result_type>(k - 1);
0111     }
0112     // Use Gaussian approximation rounded to nearest integer
0113     return result_type(sample_normal_(rng) + real_type(0.5));
0114 }
0115 //---------------------------------------------------------------------------//
0116 }  // namespace celeritas