corecel/grid/Interpolator.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/grid/Interpolator.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 #include "corecel/cont/Array.hh"
0015
0016 #include "GridTypes.hh"
0017
0018 #include "detail/InterpolatorTraits.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Interpolate, with either linear or log in x and y.
0025  *
0026  * \tparam XI Transform to apply to X coordinate
0027  * \tparam YI Transform to apply to Y coordinate
0028  * \tparam T  Floating point type
0029  *
0030  * The inputs are given as two (x,y) pairs. The same two pairs can be used to
0031  * interpolate successively with this functor.
0032  *
0033  * Interpolation on the transformed coordinates is the solution for
0034  * \f$ y \f$ in \f[
0035  \frac{f_y(y) - f_y(y_l)}{f_y(y_r) - f_y(y_l)}
0036  = \frac{f_x(x) - f_x(x_l)}{f_x(x_r) - f_x(x_l)}
0037  \f]
0038  *
0039  * where \f$ f_d(v) = v \f$ for linear interpolation on the \f$ d \f$ axis
0040  * and \f$ f_d(v) = \log v \f$ for log interpolation.
0041  *
0042  * Solving for \f$y\f$, the interpolation can be rewritten to minimize
0043  * transcendatal operations and efficiently perform multiple interpolations
0044  * over the same point range:
0045  * \f[
0046   y = f^{-1}_y \left( f_y(y_l) + \frac{ p_y(n_y(y_l), y_r) }{ p_x(n_x(x_l),
0047  x_r)} \times p_x(n_x(x_l), x) \right)
0048  \f]
0049
0050  * where transformed addition is \f$p_y(y_1, y_2) \equiv f_y(y_1) + f_y(y_2)\f$
0051  ,
0052  * transformed negation is  \f$n_y(y_1) \equiv f^{-1}_y( - f_y(y_1) )\f$
0053  * and \f$ f_y(y) = y \f$ for linear interpolation in y
0054  * and \f$ f_y(y) = \log y \f$ for log interpolation in y.
0055  *
0056  * Instantiating the interpolator precalculates the transformed intercept and
0057  * slope terms, as well as the negated x-left term.
0058  * At each evaluation of the instantiated Interpolator, only the
0059  * inverse-transform and add-transformed operation need be applied.
0060  */
0061 template<Interp XI = Interp::linear,
0062          Interp YI = Interp::linear,
0063          typename T = ::celeritas::real_type>
0064 class Interpolator
0065 {
0066   public:
0067     //!@{
0068     //! \name Type aliases
0069     using real_type = T;
0070     using Point = Array<T, 2>;
0071     //!@}
0072
0073   public:
0074     // Construct with left and right values for x and y
0075     inline CELER_FUNCTION Interpolator(Point left, Point right);
0076
0077     // Interpolate
0078     inline CELER_FUNCTION real_type operator()(real_type x) const;
0079
0080   private:
0081     real_type intercept_;  //!> f_y(y_l)
0082     real_type slope_;  //!> ratio of g(y) to g(x)
0083     real_type offset_;  //!> n_x(x_l)
0084
0085     using XTraits_t = detail::InterpolatorTraits<XI, real_type>;
0086     using YTraits_t = detail::InterpolatorTraits<YI, real_type>;
0087 };
0088
0089 //! Linear interpolation
0090 template<class T>
0091 using LinearInterpolator = Interpolator<Interp::linear, Interp::linear, T>;
0092
0093 //---------------------------------------------------------------------------//
0094 // INLINE DEFINITIONS
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Construct with left and right values for x and y.
0098  */
0099 template<Interp XI, Interp YI, class T>
0100 CELER_FUNCTION Interpolator<XI, YI, T>::Interpolator(Point left, Point right)
0101 {
0102     enum
0103     {
0104         X = 0,
0105         Y = 1
0106     };
0107
0108     CELER_EXPECT(left[X] < right[X]);
0109     CELER_EXPECT(XTraits_t::valid_domain(left[X]));
0110     CELER_EXPECT(XTraits_t::valid_domain(right[X]));
0111     CELER_EXPECT(YTraits_t::valid_domain(left[Y]));
0112     CELER_EXPECT(YTraits_t::valid_domain(right[Y]));
0113
0114     intercept_ = YTraits_t::transform(left[Y]);
0115     slope_ = (YTraits_t::add_transformed(
0116                   YTraits_t::negate_transformed(left[Y]), right[Y])
0117               / XTraits_t::add_transformed(
0118                   XTraits_t::negate_transformed(left[X]), right[X]));
0119     offset_ = XTraits_t::negate_transformed(left[X]);
0120     CELER_ENSURE(!std::isnan(intercept_) && !std::isnan(slope_)
0121                  && !std::isnan(offset_));
0122 }
0123
0124 //---------------------------------------------------------------------------//
0125 /*!
0126  * Interpolate linearly on the transformed type.
0127  */
0128 template<Interp XI, Interp YI, class T>
0129 CELER_FUNCTION auto
0130 Interpolator<XI, YI, T>::operator()(real_type x) const -> real_type
0131 {
0132     CELER_EXPECT(XTraits_t::valid_domain(x));
0133     real_type result = YTraits_t::transform_inv(
0134         std::fma(slope_, XTraits_t::add_transformed(offset_, x), intercept_));
0135
0136     CELER_ENSURE(!std::isnan(result));
0137     return result;
0138 }
0139
0140 //---------------------------------------------------------------------------//
0141 }  // namespace celeritas