random/distribution/PowerDistribution.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/PowerDistribution.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010 #include <type_traits>
0011
0012 #include "corecel/Assert.hh"
0013 #include "corecel/Macros.hh"
0014 #include "corecel/Types.hh"
0015 #include "corecel/math/Algorithms.hh"
0016
0017 #include "UniformRealDistribution.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Sample a power distribution for powers not equal to -1.
0024  *
0025  * This distribution for a power \f$ p != -1 \f$ is defined on a positive range
0026  * \f$ [a, b) \f$ and has the normalized PDF:
0027  * \f[
0028    f(x; p, a, b) = x^{p}\frac{p + 1}{b^{p + 1} - a^{p + 1}}
0029    \quad \mathrm{for } a < x < b
0030    \f]
0031  * which integrated into a CDF and inverted gives a sample:
0032  * \f[
0033     x = \left[ (b^{p + 1} - a^{p + 1})\xi + a^{p + 1} \right]^{1/(p + 1)}
0034    \f]
0035  *
0036  * The quantity in brackets is a uniform sample on
0037  * \f$ [a^{p + 1}, b^{p + 1}) \f$
0038  *
0039  * See C13A (and C15A) from \citet{everett-montecarlo-1972,
0040  * https://doi.org/10.2172/4589395} , with \f$ x \gets u \f$,
0041  * \f$ p \gets m - 1 \f$.
0042  *
0043  * \note For \f$ p = -1 \f$ see \c ReciprocalDistribution ,
0044  * and in the degenerate case of \f$ p = 0 \f$ this is mathematically
0045  * equivalent to \c UniformRealDistribution .
0046  */
0047 template<class RealType = ::celeritas::real_type>
0048 class PowerDistribution
0049 {
0050   public:
0051     //!@{
0052     //! \name Type aliases
0053     using real_type = RealType;
0054     using result_type = real_type;
0055     //!@}
0056
0057   public:
0058     // Construct on an the interval [0, 1)
0059     explicit inline CELER_FUNCTION PowerDistribution(real_type p);
0060
0061     // Construct on an arbitrary interval
0062     inline CELER_FUNCTION
0063     PowerDistribution(real_type p, real_type a, real_type b);
0064
0065     // Sample a random number according to the distribution
0066     template<class Generator>
0067     inline CELER_FUNCTION result_type operator()(Generator& rng) const;
0068
0069   private:
0070     UniformRealDistribution<real_type> sample_before_exp_;
0071     real_type exp_;
0072 };
0073
0074 //---------------------------------------------------------------------------//
0075 // INLINE DEFINITIONS
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Construct on the interval [0, 1).
0079  */
0080 template<class RealType>
0081 CELER_FUNCTION PowerDistribution<RealType>::PowerDistribution(real_type p)
0082     : sample_before_exp_{}, exp_{1 / (p + 1)}
0083 {
0084     CELER_EXPECT(p != -1);
0085 }
0086
0087 //---------------------------------------------------------------------------//
0088 /*!
0089  * Construct on the interval [a, b).
0090  *
0091  * It is allowable for the two bounds to be out of order.
0092  */
0093 template<class RealType>
0094 CELER_FUNCTION PowerDistribution<RealType>::PowerDistribution(real_type p,
0095                                                               real_type a,
0096                                                               real_type b)
0097     : sample_before_exp_{fastpow(a, p + 1), fastpow(b, p + 1)}
0098     , exp_{1 / (p + 1)}
0099 {
0100     CELER_EXPECT(p != -1);
0101     CELER_EXPECT(a >= 0);
0102     CELER_EXPECT(b >= 0);
0103 }
0104
0105 //---------------------------------------------------------------------------//
0106 /*!
0107  * Sample a random number according to the distribution.
0108  */
0109 template<class RealType>
0110 template<class Generator>
0111 CELER_FUNCTION auto
0112 PowerDistribution<RealType>::operator()(Generator& rng) const -> result_type
0113 {
0114     real_type xi = sample_before_exp_(rng);
0115     return fastpow(xi, exp_);
0116 }
0117
0118 //---------------------------------------------------------------------------//
0119 }  // namespace celeritas