Linear Interpolation for sorted arrays

Question

I have made a linear interpolation functions as a side project of mine. It assumes everything is sorted before hand - x and f(x) are the same length.

I would like to ask for:

• general recommendations
• if I can improve the design (generalize it to all iterables like std::array etc...)
• and for further performance optimizations.

#include <algorithm>
#include <cmath>
#include <type_traits>
#include <vector>

#if (__cplusplus < 202002L)
template <typename T,
          typename = std::enable_if_t<std::is_floating_point<T>::value>>
T lerp(T a, T b, T t) noexcept {
  return a + t * (b - a);
}
#else
using std::lerp;
#endif

template <typename T,
          typename = std::enable_if_t<std::is_floating_point<T>::value>>
class ListLerp {
 public:
  explicit ListLerp(const std::vector<T>& xp, const std::vector<T>& yp)
      : xp(xp), yp(yp) {
    b = this->xp.begin() + 1;
  }

  [[nodiscard]] T interp(const T x) noexcept {
    // No extrapolation
    if (x <= xp.front()) return yp.front();
    if (x >= xp.back()) return yp.back();

    // b = xp.begin() + 1;

    // find b position
    while (x > *b) b++;

    // Calculate a, a_y, b_y from b
    const auto a = b - 1;
    const auto a_y = yp.begin() + (a - xp.begin());
    const auto b_y = yp.begin() + (b - xp.begin());

    const auto t = (x - *a) / (*b - *a);

    return lerp(*a_y, *b_y, t);
  };

 private:
  typename std::vector<T>::const_iterator b;
  const std::vector<T>& xp;
  const std::vector<T>& yp;
};

template <typename T,
          typename = std::enable_if_t<std::is_floating_point<T>::value>>
std::vector<T> interp(const std::vector<T>& xp, const std::vector<T>& yp,
                      const std::vector<T>& x) {
  std::vector<T> out(x.size());

  auto lLerp = ListLerp(xp, yp);

  std::transform(x.begin(), x.end(), out.begin(),
                 [&lLerp](const auto& xi) { return lLerp.interp(xi); });

  return out;
}


///////////////
int main() {
  std::vector<double> xp = {1, 2, 3};
  std::vector<double> fp = {3, 2, 0};

  for (auto &i : interp(xp, fp, {1, 1.5, 2.72})){
    std::cout << " " << i; 
  }
  std::cout << "\n";

}

Answer 1

template <typename T,
          typename = std::enable_if_t<std::is_floating_point<T>::value>>
T lerp(T a, T b, T t) noexcept;

Looks fine. I'd prefer to see static_assert than enable_if since there's no other overloads (and if there are I'd rather fine out with a compilation error).

T lerp(T a, T b, T t) noexcept {
  return a + t * (b - a);
}

How about std::fma?

Note that GCC and Clang both do this differently... They do b * t + a * (1-t) and have special cases for t at 0 and near 1. Maybe it's better to just copy their code exactly so the compiler version doesn't change the output?

---

template <typename T,
          typename = std::enable_if_t<std::is_floating_point<T>::value>>
std::vector<T> interp(const std::vector<T>& xp, const std::vector<T>& yp,
                      const std::vector<T>& x);

There's nothing in this function that directly needs T to be a floating point type. I think it's better to remove the enable_if and let the compiler show a traceback with the "real" reason T needs to be an FP.

This function takes vectors and returns a new vector. That is a fine and very common pattern, but it's not optimally efficient. If this is a hot function, you could return an iterator/some sort of range that computes the next result when it is incremented. Boost example of this kind of thing: https://www.boost.org/doc/libs/1_66_0/libs/iterator/doc/html/iterator/specialized/transform.html.

If you want to generalize this to other iterables, then don't take std::vector... instead take a templated T const& and same for ListLerp.

---

[[nodiscard]] T interp(const T x) noexcept;

Suppose I want to discard the result? Maybe I only want to store every nth item? interp is not a const member function (one might want the mutation to happen...) and T is not an "undiscardable" type. I don't see why this needs to be here.

---

I don't think you've separated your interests very well between ListLerp::interp and the free function interp.

I think there are two good options:

• make ListLerp::interp a const member function and get rid of ListLerp::b. Then the free function interp just does a "pure" transformation of a single list.

• iterate all three lists in a single function (my preference).

E.g.

auto interp(vec const& a, vec const& b, vec const& t) {
    for idx = each index in a, b {
        ... Possibly emit a value corresponding to t ...
    }
}

This is very basic pseudo code, but the key idea is you iterate everything in one go rather than doing some iteration in a free function and some iteration in a class that holds potentially dangerous references.

---

const std::vector<T>& xp;

This (and the lines like it) is moderately dangerous since you need to ensure the vector that is passed in outlives the class.

E.g. you cannot do this:

auto fn() {
     vector<double> v;
     ListLerp ll(v); // take reference to local variable v
     return ll; // v goes out of scope; ll has dangling reference to dead v
}

Often times the best option is to take the const& and be careful about how you use it, but in this case I think it's easy enough to compute the full result in a function, not expose a class that can hold a reference to a ctor argument, and make it impossible to end up with a dangling reference.

---

It assumes everything is sorted before hand - x and f(x) are the same length.

How about asserting this?

---

// find b position
while (x > *b) b++;

If your input is sorted, this could be a binary search. May or may not be better.

---

I think your calculation at the end would be simpler if you use indices rather than iterators. Many of the lines are basically "find the iterator in A that has the same index as this iterator in B" ... easier to say "A[i]"