3
\$\begingroup\$

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

}
\$\endgroup\$
4
\$\begingroup\$
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]"

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – Stavros Avramidis May 11 at 6:35
  • \$\begingroup\$ I made a few edits in the answer. Lmk what you think. \$\endgroup\$ – sudo rm -rf slash May 11 at 12:44
  • \$\begingroup\$ I see what you mean now, thx \$\endgroup\$ – Stavros Avramidis May 11 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.