celeritas/grid/InverseCdfFinder.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 celeritas/grid/InverseCdfFinder.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/cont/Range.hh"
0012 #include "corecel/grid/Interpolator.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/SoftEqual.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Given a sampled CDF value, find the corresponding grid value.
0021  *
0022  * \tparam G Grid, e.g. \c UniformGrid or \c NonUniformGrid
0023  * \tparam C Calculate the CDF at a given grid index
0024  *
0025  * Both the input grid and the CDF must be monotonically increasing. The
0026  * sampled CDF value must be in range.
0027  */
0028 template<class G, class C>
0029 class InverseCdfFinder
0030 {
0031   public:
0032     // Construct from grid and CDF calculator
0033     inline CELER_FUNCTION InverseCdfFinder(G&& grid, C&& calc_cdf);
0034
0035     // Find and interpolate the grid value corresponding to the given CDF
0036     inline CELER_FUNCTION real_type operator()(real_type cdf) const;
0037
0038   private:
0039     G grid_;
0040     C calc_cdf_;
0041 };
0042
0043 //---------------------------------------------------------------------------//
0044 // INLINE DEFINITIONS
0045 //---------------------------------------------------------------------------//
0046 /*!
0047  * Construct from grid and CDF calculator.
0048  */
0049 template<class G, class C>
0050 CELER_FUNCTION InverseCdfFinder<G, C>::InverseCdfFinder(G&& grid, C&& calc_cdf)
0051     : grid_(celeritas::forward<G>(grid))
0052     , calc_cdf_(celeritas::forward<C>(calc_cdf))
0053 {
0054     CELER_EXPECT(grid_.size() >= 2);
0055     CELER_EXPECT(calc_cdf_[0] == 0
0056                  && soft_equal(calc_cdf_[grid_.size() - 1], real_type(1)));
0057 }
0058
0059 //---------------------------------------------------------------------------//
0060 /*!
0061  * Find and interpolate the grid value corresponding to the given CDF.
0062  */
0063 template<class G, class C>
0064 CELER_FUNCTION real_type InverseCdfFinder<G, C>::operator()(real_type cdf) const
0065 {
0066     CELER_EXPECT(cdf >= 0 && cdf < 1);
0067
0068     // Find the grid index of the sampled CDF value
0069     Range indices(grid_.size() - 1);
0070     auto iter = celeritas::lower_bound(
0071         indices.begin(), indices.end(), cdf, [this](size_type i, real_type c) {
0072             return calc_cdf_[i] < c;
0073         });
0074     if (calc_cdf_[*iter] != cdf)
0075     {
0076         --iter;
0077     }
0078     size_type i = iter - indices.begin();
0079
0080     // Calculate the grid value corresponding to the sampled CDF value
0081     return LinearInterpolator<real_type>{
0082         {calc_cdf_[i], grid_[i]}, {calc_cdf_[i + 1], grid_[i + 1]}}(cdf);
0083 }
0084
0085 //---------------------------------------------------------------------------//
0086 }  // namespace celeritas