Linear Interpolation for sorted arrays

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";

}


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]"

• Thanks for your answer. I will play around with the suggestions. I have 2 questions. One about your note in the separation between ListLerp::interp - if you could give me an example. The reason ListLerp exists as a class in general is not to do a linear search in every item. Also if you could explain the co st reference issue a bit. Commented May 11, 2020 at 6:35
• I made a few edits in the answer. Lmk what you think. Commented May 11, 2020 at 12:44
• I see what you mean now, thx Commented May 11, 2020 at 14:35