Interpolator.hpp

#pragma once
#include <cassert>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_version.h>
#include <iostream>
#include <map>
#include <vector>
 
// This should make code work with old versions of GSL
#ifdef GSL_MAJOR_VERSION
#if GSL_MAJOR_VERSION == 1
#define GSL_VERSION_1
#endif
#endif
 
namespace Interpolator {
 
 
enum class Method {
  linear,
  polynomial,
  cspline,
  cspline_periodic,
  akima,
  akima_periodic,
  steffen
};
 
static const std::map<Method, const gsl_interp_type *> interp_method{
    {Method::linear, gsl_interp_linear},
    {Method::polynomial, gsl_interp_polynomial},
    {Method::cspline, gsl_interp_cspline},
    {Method::cspline_periodic, gsl_interp_cspline_periodic},
    {Method::akima, gsl_interp_akima},
    {Method::akima_periodic, gsl_interp_akima_periodic}
#ifndef GSL_VERSION_1
    ,
    {Method::steffen, gsl_interp_steffen}
#endif
};
 
//==============================================================================
 
 
class Interp {
 
private:
  gsl_interp_accel *acc;
  gsl_spline *spline;
  double x0, xmax;
 
public:
  Interp(const std::vector<double> &x, const std::vector<double> &y,
         Method method = Method::cspline)
      : acc(gsl_interp_accel_alloc()),
        spline(gsl_spline_alloc(interp_method.at(method), x.size())),
        x0(x.front()),
        xmax(x.back()) {
    assert(x.size() == y.size() &&
           "In Interp, input vectors x and y must have same size");
    assert(x.size() >= gsl_interp_type_min_size(interp_method.at(method)) &&
           "In Interp, certain interpolation methods require a minimum number "
           "of points");
    gsl_spline_init(spline, x.data(), y.data(), x.size());
  }
  ~Interp() {
    gsl_spline_free(spline);
    gsl_interp_accel_free(acc);
  }
  Interp(const Interpolator::Interp &) = delete;
  Interp &operator=(const Interpolator::Interp &) = delete;
 
  double interp(double x) const {
    if (x < x0 || x > xmax)
      return 0.0;
    else
      return gsl_spline_eval(spline, x, acc);
  }
 
  std::vector<double> interp(const std::vector<double> &x) const {
    std::vector<double> y;
    y.reserve(x.size());
    for (auto xi : x)
      y.push_back(interp(xi));
    return y;
  }
 
  double operator()(double x) const { return interp(x); }
 
  std::vector<double> operator()(const std::vector<double> &x) const {
    return interp(x);
  }
};
 
//==============================================================================
 
 
inline std::vector<double> interpolate(const std::vector<double> &x_in,
                                       const std::vector<double> &y_in,
                                       const std::vector<double> &x_out,
                                       Method method = Method::cspline) {
  Interp i_func(x_in, y_in, method);
  return i_func(x_out);
}
 
//==============================================================================
 
static constexpr bool has_steffen_method() {
#ifndef GSL_VERSION_1
  return true;
#else
  return false;
#endif
}
 
} // namespace Interpolator