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