Interpolate.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_
#include <molpro/linalg/scalar_traits.h>
 
#include <cmath>
#include <limits>
#include <ostream>
#include <string>
#include <vector>
 
namespace molpro::linalg::itsolv {
 
template <typename value_type = double>
class Interpolator {
  static_assert(!molpro::linalg::is_complex<value_type>{},
                "the interpolant orders function values in order to minimise them, so it is defined "
                "only for a real type; pass the real type underlying your scalar, as Optimize does");
 
public:
  using value_t = value_type;
 
  struct point {
    value_type x;                                                     //< abscissa
    value_type f = std::numeric_limits<value_type>::quiet_NaN();      //< function value at x
    value_type f1 = std::numeric_limits<value_type>::quiet_NaN();     //< function first gradient at x
    value_type f2 = std::numeric_limits<value_type>::quiet_NaN();     //< function second gradient at x
  };
  explicit Interpolator(point p0, point p1, std::string interpolant = "cubic", int verbosity = 0);
  point operator()(value_type x) const;
  point minimize(value_type xa, value_type xb, size_t bracket_grid = 100, size_t max_bracket_grid = 100000,
                 bool analytic = true) const;
  point minimize_cubic() const;
  static std::vector<std::string> interpolants();
 
  const std::vector<value_type>& parameters() const { return m_parameters; }
 
private:
  const point m_p0, m_p1;
  const std::string m_interpolant;
  std::vector<value_type> m_parameters;
};
 
using Interpolate = Interpolator<double>;
 
template <typename value_type>
bool operator==(const typename Interpolator<value_type>::point& lhs,
                const typename Interpolator<value_type>::point& rhs) {
  return lhs.x == rhs.x && lhs.f == rhs.f && lhs.f1 == rhs.f1;
}
 
inline bool operator==(const Interpolate::point& lhs, const Interpolate::point& rhs) {
  return lhs.x == rhs.x && lhs.f == rhs.f && lhs.f1 == rhs.f1;
}
 
template <typename value_type>
std::ostream& operator<<(std::ostream& os, const Interpolator<value_type>& interpolant) {
  for (const auto& parameter : interpolant.parameters())
    os << " " << parameter;
  return os;
}
 
template <typename value_type>
std::ostream& operator<<(std::ostream& os, const typename Interpolator<value_type>::point& p) {
  os << "x=" << p.x << ", value=" << p.f << ", gradient=" << p.f1 << ", curvature=" << p.f2;
  return os;
}
 
inline std::ostream& operator<<(std::ostream& os, const Interpolate::point& p) {
  os << "x=" << p.x << ", value=" << p.f << ", gradient=" << p.f1 << ", curvature=" << p.f2;
  return os;
}
 
extern template class Interpolator<double>;
extern template class Interpolator<long double>;
 
} // namespace molpro::linalg::itsolv
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_