corecel/grid/Interpolator.hh

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