cubic_spline.h

Go to the documentation of this file.

#ifndef CUBIC_SPLINE_INTERPOLATION_H_
#define CUBIC_SPLINE_INTERPOLATION_H_
 
#include <vector>
#include <cstddef>
 
namespace ModuleBase
{
 
class CubicSpline
{    
    //*****************************************************************
    //                      boundary condition
    //*****************************************************************
 
public:
 
    enum class BoundaryType
    {
        not_a_knot,
        first_deriv,
        second_deriv,
        periodic
    };
 
 
    struct BoundaryCondition
    {
        // for not_a_knot and periodic
        BoundaryCondition(BoundaryType type = BoundaryType::not_a_knot);
 
        // for first/second_deriv
        BoundaryCondition(BoundaryType type, double val);
 
        BoundaryType type;
        double val = 0.0;
    };
 
 
    //*****************************************************************
    //                      interpolant object
    //*****************************************************************
 
public:
 
    CubicSpline()                   = delete;
    CubicSpline(CubicSpline const&) = default;
    CubicSpline(CubicSpline &&)     = default;
 
    CubicSpline& operator=(CubicSpline const&)  = default;
    CubicSpline& operator=(CubicSpline &&)      = default;
 
    ~CubicSpline() = default; 
 
 
    CubicSpline(
        int n,
        const double* x,
        const double* y,
        const BoundaryCondition& bc_start = {},
        const BoundaryCondition& bc_end = {}
    );
 
 
    CubicSpline(
        int n,
        double x0,
        double dx,
        const double* y,
        const BoundaryCondition& bc_start = {},
        const BoundaryCondition& bc_end = {}
    );
 
 
    CubicSpline(int n, const double* x);
 
    CubicSpline(int n, double x0, double dx);
 
    void add(
        const double* y,
        const BoundaryCondition& bc_start = {},
        const BoundaryCondition& bc_end = {}
    );
 
 
    void eval(
        int n_interp,
        const double* x_interp,
        double* y_interp,
        double* dy_interp = nullptr,
        double* d2y_interp = nullptr,
        int i_spline = 0
    ) const;
 
 
    void multi_eval(
        int n_spline,
        const int* i_spline,
        double x_interp,
        double* y_interp,
        double* dy_interp = nullptr,
        double* d2y_interp = nullptr
    ) const;
 
 
    void multi_eval(
        double x_interp,
        double* y_interp,
        double* dy_interp = nullptr,
        double* d2y_interp = nullptr
    ) const;
 
 
    void reserve(int n_spline) { y_.reserve(n_spline * n_ * 2); }
 
 
    size_t heap_usage() const { return (x_.capacity() + y_.capacity()) * sizeof(double); }
 
    double xmin() const { return xmin_; }
 
    double xmax() const { return xmax_; }
 
    int n_spline() const { return n_spline_; }
 
 
private:
 
    int n_spline_ = 0;
 
    int n_ = 0;
 
    double xmin_ = 0.0;
 
    double xmax_ = 0.0;
 
    double dx_ = 0.0;
 
    std::vector<double> x_;
 
    std::vector<double> y_;
 
 
    //*****************************************************************
    //                      static functions
    //*****************************************************************
 
public:
 
    static void build(
        int n,
        const double* x,
        const double* y,
        const BoundaryCondition& bc_start,
        const BoundaryCondition& bc_end,
        double* dy
    );
 
 
    static void build(
        int n,
        double dx,
        const double* y,
        const BoundaryCondition& bc_start,
        const BoundaryCondition& bc_end,
        double* dy
    );
 
 
    static void eval(
        int n,
        const double* x,
        const double* y,
        const double* dy,
        int n_interp,
        const double* x_interp,
        double* y_interp,
        double* dy_interp = nullptr,
        double* d2y_interp = nullptr
    );
 
 
    static void eval(
        int n,
        double x0,
        double dx,
        const double* y,
        const double* dy,
        int n_interp,
        const double* x_interp,
        double* y_interp,
        double* dy_interp = nullptr,
        double* d2y_interp = nullptr
    );
 
 
private:
 
    static void _build(
        int n,
        const double* dx,
        const double* y,
        const BoundaryCondition& bc_start,
        const BoundaryCondition& bc_end,
        double* dy
    );
 
 
    static inline int _index(int n, const double* x, double target);
 
 
    static inline int _index(int n, double x0, double dx, double target);
 
 
    static inline void _cubic(
        int n,
        const double* w,
        const double* c0,
        const double* c1,
        const double* c2,
        const double* c3,
        double* y,
        double* dy,
        double* d2y
    );
 
 
    static void _validate_build(
        int n,
        const double* dx,
        const double* y,
        const BoundaryCondition& bc_start,
        const BoundaryCondition& bc_end
    );
 
 
    static void _validate_eval(
        int n,
        const double (&u)[2],
        const double* x,
        const double* y,
        const double* dy,
        int n_interp,
        const double* x_interp
    );
 
 
    static void _solve_cyctri(int n, double* d, double* u, double* l, double* b);
};
 
} // namespace ModuleBase
 
#endif