.. _program_listing_file_src_interpolation_Interpolator1D.hpp:

Program Listing for File Interpolator1D.hpp
===========================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_interpolation_Interpolator1D.hpp>` (``src/interpolation/Interpolator1D.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /*
    * @file Interpolator1D.hpp
    * @author Mark Krumholz
    * @brief Provides a class to represent a single 1D interpolator
    * @date 2024-06-27
    * @copyright Copyright (c) 2026 Mark Krumholz. All rights reserved.
   */
   
   #ifndef INTERPOLATOR1D_HPP
   #define INTERPOLATOR1D_HPP
   
   #include "../utils/ThreadVec.hpp"
   #include <algorithm>
   #include <array>
   #include <cassert>
   #include <cstddef>
   #include <gsl/gsl_interp.h>
   #include <ranges>
   #include <sstream>
   #include <stdexcept>
   #include <utility>
   #include <vector>
   
   // Disable linting for array bounds checking in this
   // file, since the overhead associated with enforcing
   // such checks severely interferes with performance
   // NOLINTBEGIN(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
   
   namespace interp
   {
   
       template <size_t NF = 1>
       class Interpolator1D
       {
       public:
   
           // Constructors and destructor
   
           Interpolator1D(
               const std::vector<double>& x,
               const std::array<std::vector<double>, NF>& f,
               const gsl_interp_type* interpType = gsl_interp_steffen
           ) :
           x_(x),
           f_(f),
           interp_({nullptr})
           {
               checkSorted();
               cleanDuplicates();
               checkMinSize(interpType);
               interpInit(interpType);
           }
   
           Interpolator1D(
               const std::vector<double>& x,
               const std::vector<double>& f,
               const gsl_interp_type* interpType = gsl_interp_steffen
           ) requires (NF == 1)
           : 
           Interpolator1D(x, 
               std::array<std::vector<double>, 1>({f}),
               interpType)
           { }
   
           // NOLINTNEXTLINE(portability-template-virtual-member-function) -- only reachable through unique_ptr<Interpolator1D<NF>> (e.g. specsyn::Specsyn::Isochrone), never a raw new/delete pair the compiler could implicitly instantiate differently; this project also always builds a given binary with a single, consistent compiler, so cross-compiler ABI divergence in the implicit-instantiation rules for this destructor is not a real risk here
           virtual ~Interpolator1D()
           {
               // Free gsl data
               for (auto& i : interp_) { gsl_interp_free(i); }
               for (auto& a : acc_) { gsl_interp_accel_free(a); }
           }
   
           // Disallow move and copy constructors, since these would make
           // copies of gsl opaque objects that we can't manage with
           // smart pointers, resulting in potential memory corruption
           // when two Interpolator1D objects pointing to the same gsl
           // opaque objects pass out of scope and try to free the same
           // memory. We manage this by only ever dealing with unique_ptr's
           // to Interpolator1D objects elsewhere in the code, to ensure
           // that the destructor is invoked once and only once.
           Interpolator1D(const Interpolator1D&) = delete;
           auto operator=(const Interpolator1D&) -> Interpolator1D& = delete;
           Interpolator1D(Interpolator1D&&) = delete;
           auto operator=(Interpolator1D&&) -> Interpolator1D& = delete;
   
           [[nodiscard]] auto xMin() const { return x_.front(); }
   
           [[nodiscard]] auto xMax() const { return x_.back(); }
   
           [[nodiscard]] auto xRange() const { return std::make_pair(xMin(), xMax()); }
   
           using PostProcess = void (*)(std::array<double, NF>&);
   
           void setPostProcess(const PostProcess postProcess) { postProcess_ = postProcess; }
   
           [[nodiscard]] auto operator()(double x) const
           {
               assert(x >= x_.front() && x <= x_.back());
               std::array<double, NF> result = {};
               for (size_t i = 0; i < NF; ++i)
               {
                   result[i] = gsl_interp_eval(interp_[i], x_.data(), // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i is a loop index bounded by compile-time constant NF
                       f_[i].data(), x, acc_()); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i is a loop index bounded by compile-time constant NF
               }
               if (postProcess_ != nullptr) { postProcess_(result); }
               if constexpr(NF > 1) { return result; }
               else { return result[0]; }
           }
   
           [[nodiscard]] auto operator()(double x, size_t idx) const
           {
               assert(x >= x_.front() && x <= x_.back());
               assert(idx < NF);
               if (postProcess_ != nullptr)
               {
                   if constexpr(NF > 1) { return (*this)(x)[idx]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- idx is asserted < NF just above
                   else { return (*this)(x); }
               }
               return evalRaw(x, idx);
           }
   
           [[nodiscard]] auto integ(double x0, double x1) const
           {
               assert(x0 >= x_.front() && x0 <= x_.back());
               assert(x1 >= x_.front() && x1 <= x_.back());
               assert(x0 <= x1);
               std::array<double, NF> result = {};
               for (size_t i = 0; i < NF; ++i)
               {
                   result[i] = gsl_interp_eval_integ(interp_[i], x_.data(), // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i is a loop index bounded by compile-time constant NF
                       f_[i].data(), x0, x1, acc_()); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i is a loop index bounded by compile-time constant NF
               }
               if constexpr(NF > 1) { return result; }
               else { return result[0]; }
           }
   
           [[nodiscard]] auto integ(double x0, double x1, size_t idx) const
           {
               assert(x0 >= x_.front() && x0 <= x_.back());
               assert(x1 >= x_.front() && x1 <= x_.back());
               assert(x0 <= x1);
               assert(idx < NF);
               return gsl_interp_eval_integ(interp_[idx], x_.data(), f_[idx].data(), x0, x1, acc_()); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- idx is asserted < NF just above
           }
   
           [[nodiscard]] auto integ(
               const Interpolator1D& g, double x0, double x1, size_t idx) const -> double
           {
               // NF equality is enforced at compile time: both are Interpolator1D<NF>
               assert(idx < NF);
   
               const double xLo = std::max({x_.front(), g.x_.front(), x0});
               const double xHi = std::min({x_.back(),  g.x_.back(),  x1});
               if (xLo >= xHi) { return 0.0; }
   
               // Merged, deduplicated x-grid: xLo, all internal points from either
               // grid that fall strictly inside (xLo, xHi), then xHi
               std::vector<double> xGrid;
               xGrid.reserve(x_.size() + g.x_.size() + 2);
               xGrid.push_back(xLo);
               for (const double xi : x_)   { if (xi > xLo && xi < xHi) { xGrid.push_back(xi); } }
               for (const double xi : g.x_) { if (xi > xLo && xi < xHi) { xGrid.push_back(xi); } }
               xGrid.push_back(xHi);
               std::ranges::sort(xGrid);
               auto toErase = std::ranges::unique(xGrid);
               xGrid.erase(toErase.begin(), toErase.end());
   
               // Evaluate f_idx(x) * g_idx(x) at each grid point
               std::vector<double> fgGrid(xGrid.size());
               for (size_t j = 0; j < xGrid.size(); ++j)
               {
                   fgGrid[j] = evalRaw(xGrid[j], idx) * g.evalRaw(xGrid[j], idx);
               }
   
               // Integrate the product interpolant
               const Interpolator1D<1> fgInterp(xGrid, fgGrid, gsl_interp_linear);
               return fgInterp.integ(xLo, xHi);
           }
   
           [[nodiscard]] auto integ(
               const Interpolator1D& g, double x0, double x1) const
           {
               // NF equality is enforced at compile time: both are Interpolator1D<NF>
               const double xLo = std::max({x_.front(), g.x_.front(), x0});
               const double xHi = std::min({x_.back(),  g.x_.back(),  x1});
               if (xLo >= xHi)
               {
                   if constexpr (NF > 1) { return std::array<double, NF>{}; }
                   else { return 0.0; }
               }
   
               // Merged, deduplicated x-grid: xLo, all internal points from either
               // grid that fall strictly inside (xLo, xHi), then xHi
               std::vector<double> xGrid;
               xGrid.reserve(x_.size() + g.x_.size() + 2);
               xGrid.push_back(xLo);
               for (const double xi : x_)   { if (xi > xLo && xi < xHi) { xGrid.push_back(xi); } }
               for (const double xi : g.x_) { if (xi > xLo && xi < xHi) { xGrid.push_back(xi); } }
               xGrid.push_back(xHi);
               std::ranges::sort(xGrid);
               auto toErase = std::ranges::unique(xGrid);
               xGrid.erase(toErase.begin(), toErase.end());
   
               // Evaluate f_i(x) * g_i(x) at each grid point for every field
               std::array<std::vector<double>, NF> fgGrid;
               for (auto& v : fgGrid) { v.resize(xGrid.size()); }
               for (size_t j = 0; j < xGrid.size(); ++j)
               {
                   for (size_t i = 0; i < NF; ++i)
                   {
                       fgGrid[i][j] = evalRaw(xGrid[j], i) * g.evalRaw(xGrid[j], i); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i is bounded by compile-time constant NF
                   }
               }
   
               // Integrate the product interpolant
               const Interpolator1D<NF> fgInterp(xGrid, fgGrid, gsl_interp_linear);
               return fgInterp.integ(xLo, xHi);
           }
   
       private:
   
           void checkSorted()
           {
               static_assert(NF > 0, "Interpolator1D: number of quantities to interpolate must be >= 1");
               if (!std::ranges::is_sorted(x_))
               {
                   throw std::runtime_error("Interpolator1D: x must be strictly increasing");
               }
           }
   
           void checkMinSize(const gsl_interp_type *interpType)
           {
               if (x_.size() < gsl_interp_type_min_size(interpType))
               {
                   std::stringstream ss;
                   ss << "Interpolator1D: interpolation type requires minimum size "
                       << gsl_interp_type_min_size(interpType)
                       << ", but only " << x_.size() << " points provided";
                   throw std::runtime_error(ss.str());
               }
           }
   
           void cleanDuplicates()
           {
               if (x_.size() < 2) { return; }
               size_t w = 0;
               for (size_t r = 1; r < x_.size(); ++r)
               {
                   if (x_[r] == x_[w]) { continue; }
                   ++w;
                   x_[w] = x_[r];
                   for (auto &fi : f_) { fi[w] = fi[r]; }
               }
               const size_t newSize = w + 1;
               x_.resize(newSize);
               for (auto &fi : f_) { fi.resize(newSize); }
           }
   
           [[nodiscard]] auto evalRaw(const double x, const size_t idx) const -> double
           {
               assert(idx < NF);
               return gsl_interp_eval(interp_[idx], x_.data(), f_[idx].data(), x, acc_()); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- idx is asserted < NF just above
           }
   
           void interpInit(const gsl_interp_type *interpType)
           {
               for (auto&& [i, d] : std::views::zip(interp_, f_))
               {
                   if (x_.size() != d.size())
                   {
                       throw std::runtime_error("Interpolator1D: x and f must be of same length");
                   }
                   if (x_.size() > gsl_interp_type_min_size(interpType))
                   {
                       i = gsl_interp_alloc(interpType, x_.size());
                   }
                   else
                   {
                       i = gsl_interp_alloc(gsl_interp_linear, x_.size());
                   }
                   gsl_interp_init(i, x_.data(), d.data(), x_.size());
               }
               for (auto& a : acc_) { a = gsl_interp_accel_alloc(); }
           }
   
           // Internal storage
           std::vector<double> x_;                 
           std::array<std::vector<double>, NF> f_; 
           std::array<gsl_interp *, NF> interp_;   
           utils::ThreadVec<gsl_interp_accel *> acc_; 
           PostProcess postProcess_ = nullptr; 
       };
   
   } // namespace interp
   
   // NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
   
   #endif // INTERPOLATOR1D_HPP