Program Listing for File Interpolator1D.hpp
↰ Return to documentation for file (src/interpolation/Interpolator1D.hpp)
/*
* @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