/*
* spline.h
*
* simple cubic spline interpolation library without external
* dependencies
*
* ---------------------------------------------------------------------
* Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com)
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
* ---------------------------------------------------------------------
*
*/
#ifndef TK_SPLINE_H
#define TK_SPLINE_H
#include
#include
#include
#include
namespace tk
{
// band matrix solver
class band_matrix
{
private:
std::vector< std::vector > m_upper; // upper band
std::vector< std::vector > m_lower; // lower band
public:
inline band_matrix() = default; // constructor
inline band_matrix(int dim, int n_u, int n_l); // constructor
inline void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l
inline int dim() const; // matrix dimension
inline int num_upper() const
{
return m_upper.size()-1;
}
inline int num_lower() const
{
return m_lower.size()-1;
}
// access operator
inline double & operator () (int i, int j); // write
inline double operator () (int i, int j) const; // read
// we can store an additional diogonal (in m_lower)
inline double& saved_diag(int i);
inline double saved_diag(int i) const;
inline void lu_decompose();
inline std::vector r_solve(const std::vector& b) const;
inline std::vector l_solve(const std::vector& b) const;
inline std::vector lu_solve(const std::vector& b,
bool is_lu_decomposed=false);
};
enum class spline_type {
linear = 1,
cubic = 3,
dflt = cubic
};
enum class extrapolation_method {
linear = 1,
quadratic = 2,
dflt = quadratic
};
// spline interpolation
class spline
{
public:
enum bd_type {
first_deriv = 1,
second_deriv = 2
};
private:
// find the closest point m_x[idx] < x, idx=0 even if x m_x,m_y; // x,y coordinates of points
// interpolation parameters
// f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
std::vector m_a,m_b,m_c; // spline coefficients
double m_b0, m_c0; // for left extrapol
public:
inline spline() = default;
inline spline(std::vector xs, std::vector ys,
spline_type type = spline_type::dflt,
extrapolation_method extrapolation = extrapolation_method::dflt,
bd_type left = second_deriv, double left_value = 0.0,
bd_type right = second_deriv, double right_value = 0.0);
// optional, but if called it has to come be before set_points()
inline double operator() (double x) const;
inline double deriv(int order, double x) const;
};
// ---------------------------------------------------------------------
// implementation part, which could be separated into a cpp file
// ---------------------------------------------------------------------
// band_matrix implementation
// -------------------------
band_matrix::band_matrix(int dim, int n_u, int n_l)
{
resize(dim, n_u, n_l);
}
void band_matrix::resize(int dim, int n_u, int n_l)
{
assert(dim>0);
assert(n_u>=0);
assert(n_l>=0);
m_upper.resize(n_u+1);
m_lower.resize(n_l+1);
for(size_t i=0; i0) {
return m_upper[0].size();
} else {
return 0;
}
}
// defines the new operator (), so that we can access the elements
// by A(i,j), index going from i=0,...,dim()-1
double & band_matrix::operator () (int i, int j)
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i=0) && (j diogonal, k<0 lower left part, k>0 upper right part
if(k>=0) return m_upper[k][i];
else return m_lower[-k][i];
}
double band_matrix::operator () (int i, int j) const
{
int k=j-i; // what band is the entry
assert( (i>=0) && (i=0) && (j diogonal, k<0 lower left part, k>0 upper right part
if(k>=0) return m_upper[k][i];
else return m_lower[-k][i];
}
// second diag (used in LU decomposition), saved in m_lower
double band_matrix::saved_diag(int i) const
{
assert( (i>=0) && (i=0) && (idim(); i++) {
assert(this->operator()(i,i)!=0.0);
this->saved_diag(i)=1.0/this->operator()(i,i);
j_min=std::max(0,i-this->num_lower());
j_max=std::min(this->dim()-1,i+this->num_upper());
for(int j=j_min; j<=j_max; j++) {
this->operator()(i,j) *= this->saved_diag(i);
}
this->operator()(i,i)=1.0; // prevents rounding errors
}
// Gauss LR-Decomposition
for(int k=0; kdim(); k++) {
i_max=std::min(this->dim()-1,k+this->num_lower()); // num_lower not a mistake!
for(int i=k+1; i<=i_max; i++) {
assert(this->operator()(k,k)!=0.0);
x=-this->operator()(i,k)/this->operator()(k,k);
this->operator()(i,k)=-x; // assembly part of L
j_max=std::min(this->dim()-1,k+this->num_upper());
for(int j=k+1; j<=j_max; j++) {
// assembly part of R
this->operator()(i,j)=this->operator()(i,j)+x*this->operator()(k,j);
}
}
}
}
// solves Ly=b
std::vector band_matrix::l_solve(const std::vector& b) const
{
assert( this->dim()==(int)b.size() );
std::vector x(this->dim());
int j_start;
double sum;
for(int i=0; idim(); i++) {
sum=0;
j_start=std::max(0,i-this->num_lower());
for(int j=j_start; joperator()(i,j)*x[j];
x[i]=(b[i]*this->saved_diag(i)) - sum;
}
return x;
}
// solves Rx=y
std::vector band_matrix::r_solve(const std::vector& b) const
{
assert( this->dim()==(int)b.size() );
std::vector x(this->dim());
int j_stop;
double sum;
for(int i=this->dim()-1; i>=0; i--) {
sum=0;
j_stop=std::min(this->dim()-1,i+this->num_upper());
for(int j=i+1; j<=j_stop; j++) sum += this->operator()(i,j)*x[j];
x[i]=( b[i] - sum ) / this->operator()(i,i);
}
return x;
}
std::vector band_matrix::lu_solve(const std::vector& b,
bool is_lu_decomposed)
{
assert( this->dim()==(int)b.size() );
std::vector x,y;
if(is_lu_decomposed==false) {
this->lu_decompose();
}
y=this->l_solve(b);
x=this->r_solve(y);
return x;
}
// spline implementation
// -----------------------
inline spline::spline(std::vector xs, std::vector ys,
spline_type type,
extrapolation_method extrapolation,
bd_type left, double left_value,
bd_type right, double right_value)
: m_x(std::move(xs)), m_y(std::move(ys)) {
int n=m_x.size();
// TODO: maybe sort x and y, rather than returning an error
for(int i=0; i rhs(n);
for(int i=1; im_x[n-1]) {
// extrapolation to the right
interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1];
} else {
// interpolation
interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
}
return interpol;
}
double spline::deriv(int order, double x) const
{
assert(order>0);
const size_t n=m_x.size();
const std::size_t idx = closest_idx_to(x);
const double h=x-m_x[idx];
double interpol;
if(xm_x[n-1]) {
// extrapolation to the right
switch(order) {
case 1:
interpol=2.0*m_b[n-1]*h + m_c[n-1];
break;
case 2:
interpol=2.0*m_b[n-1];
break;
default:
interpol=0.0;
break;
}
} else {
// interpolation
switch(order) {
case 1:
interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx];
break;
case 2:
interpol=6.0*m_a[idx]*h + 2.0*m_b[idx];
break;
case 3:
interpol=6.0*m_a[idx];
break;
default:
interpol=0.0;
break;
}
}
return interpol;
}
} // namespace tk
#endif /* TK_SPLINE_H */