random/distribution/ReciprocalDistribution.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/ReciprocalDistribution.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
0015 #include "GenerateCanonical.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Reciprocal or log-uniform distribution.
0022  *
0023  * This distribution is defined on a positive range \f$ [a, b) \f$ and has the
0024  * normalized PDF:
0025  * \f[
0026    f(x; a, b) = \frac{1}{x (\ln b - \ln a)} \quad \mathrm{for} \ a \le x < b
0027    \f]
0028  * which integrated into a CDF and inverted gives a sample:
0029  * \f[
0030   x = a \left( \frac{b}{a} \right)^{\xi}
0031     = a \exp\!\left(\xi \log \frac{b}{a} \right)
0032    \f]
0033  */
0034 template<class RealType = ::celeritas::real_type>
0035 class ReciprocalDistribution
0036 {
0037   public:
0038     //!@{
0039     //! \name Type aliases
0040     using real_type = RealType;
0041     using result_type = real_type;
0042     //!@}
0043
0044   public:
0045     // Construct on an the interval [a, 1]
0046     explicit inline CELER_FUNCTION ReciprocalDistribution(real_type a);
0047
0048     // Construct on an arbitrary interval
0049     inline CELER_FUNCTION ReciprocalDistribution(real_type a, real_type b);
0050
0051     // Sample a random number according to the distribution
0052     template<class Generator>
0053     inline CELER_FUNCTION result_type operator()(Generator& rng) const;
0054
0055   private:
0056     RealType a_;
0057     RealType logratio_;
0058 };
0059
0060 //---------------------------------------------------------------------------//
0061 // INLINE DEFINITIONS
0062 //---------------------------------------------------------------------------//
0063 /*!
0064  * Construct on the interval [a, 1).
0065  *
0066  * The distribution is equivalent to switching a and b, and using
0067  * \f$ \xi' = 1 - \xi \f$.
0068  */
0069 template<class RealType>
0070 CELER_FUNCTION
0071 ReciprocalDistribution<RealType>::ReciprocalDistribution(real_type a)
0072     : ReciprocalDistribution(1, a)
0073 {
0074 }
0075
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Construct on the interval [a, b).
0079  *
0080  * It is allowable for the two bounds to be out of order.
0081  *
0082  * Note that writing as \code (1/a) * b \endcode allows the compiler to
0083  * optimize better for the constexpr case a=1.
0084  */
0085 template<class RealType>
0086 CELER_FUNCTION
0087 ReciprocalDistribution<RealType>::ReciprocalDistribution(real_type a,
0088                                                          real_type b)
0089     : a_(a), logratio_(std::log((1 / a) * b))
0090 {
0091     CELER_EXPECT(a > 0);
0092     CELER_EXPECT(b > 0);
0093 }
0094
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Sample a random number according to the distribution.
0098  */
0099 template<class RealType>
0100 template<class Generator>
0101 CELER_FUNCTION auto
0102 ReciprocalDistribution<RealType>::operator()(Generator& rng) const -> result_type
0103 {
0104     return a_ * std::exp(logratio_ * generate_canonical<RealType>(rng));
0105 }
0106
0107 //---------------------------------------------------------------------------//
0108 }  // namespace celeritas