Interpolate.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_
#include <cmath>
#include <string>
#include <vector>
 
namespace molpro::linalg::itsolv {
 
class Interpolate {
public:
  struct point {
    double x;                      //< abscissa
    double f = std::nan("unset");  //< function value at x
    double f1 = std::nan("unset"); //< function first gradient at x
    double f2 = std::nan("unset"); //< function second gradient at x
    friend std::ostream& operator<<(std::ostream& os, const point& p);
  };
  explicit Interpolate(point p0, point p1, std::string interpolant = "cubic", int verbosity = 0);
  point operator()(double x) const;
  Interpolate::point minimize(double xa, double xb, size_t bracket_grid = 100, size_t max_bracket_grid = 100000, bool analytic=true) const;
  Interpolate::point minimize_cubic() const;
  static std::vector<std::string> interpolants();
 
  const std::vector<double>& parameters() const { return m_parameters; }
  friend std::ostream& operator<<(std::ostream& os, const Interpolate& interpolant);
 
private:
  const point m_p0, m_p1;
  const std::string m_interpolant;
  double m_c0, m_c1, m_c2, m_c3; //< parameters defining the interpolant
  std::vector<double> m_parameters;
};
 
inline bool operator==(const Interpolate::point& lhs, const Interpolate::point& rhs) {
  return lhs.x == rhs.x && lhs.f == rhs.f && lhs.f1 == rhs.f1;
}
std::ostream& operator<<(std::ostream& os, const Interpolate& interpolant);
std::ostream& operator<<(std::ostream& os, const Interpolate::point& p);
 
} // namespace molpro::linalg::itsolv
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_INTERPOLATE_H_