random/distribution/PoissonDistribution.hh

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