89 lines
1.8 KiB
C++
89 lines
1.8 KiB
C++
#pragma once
|
|
|
|
#include "./numerics/interpolation1d/interpolation1d_base.h"
|
|
|
|
#include "./utils/vector.h"
|
|
#include "./numerics/min.h"
|
|
#include "./numerics/max.h"
|
|
|
|
|
|
namespace numerics{
|
|
|
|
template <typename T>
|
|
struct interp_barycentric : Base_interp<T> {
|
|
using Base = Base_interp<T>;
|
|
|
|
// bring base data members into scope (or use this->xx / this->yy below)
|
|
using Base::xx;
|
|
using Base::yy;
|
|
using Base::n;
|
|
|
|
utils::Vector<T> w;
|
|
int64_t d;
|
|
|
|
interp_barycentric(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t dd)
|
|
: Base_interp<T>(xv, &yv[0], xv.size()), w(n,T{0}), d(dd) {
|
|
// Constructor arguments are x and y vectors of length n, and order d of desired approximation.
|
|
if (n <= d){
|
|
throw std::invalid_argument("d too large for number of points in interp_barycentric");
|
|
}
|
|
for (int64_t k = 0; k < n; ++k){
|
|
int64_t imin = numerics::max(k-d, static_cast<int64_t>(0));
|
|
int64_t imax;
|
|
if (k >= n - d) {
|
|
imax = n - d - 1;
|
|
} else {
|
|
imax = k;
|
|
}
|
|
T temp;
|
|
if ( (imin & 1) != 0 ) { // odd?
|
|
temp = T{-1};
|
|
} else { // even
|
|
temp = T{1};
|
|
}
|
|
T sum = T{0};
|
|
|
|
for (int64_t i = imin; i <= imax; ++i){
|
|
int64_t jmax = numerics::min(i+d, n-1);
|
|
T term = T{1};
|
|
for (int64_t j = i; j <= jmax; ++j){
|
|
if (j == k){
|
|
continue;
|
|
}
|
|
term *= (xx[k] - xx[j]);
|
|
}
|
|
term = temp/term;
|
|
temp = -temp;
|
|
sum += term;
|
|
}
|
|
w[k] = sum;
|
|
}
|
|
}
|
|
|
|
T rawinterp(int64_t jl, T x) override{
|
|
|
|
T num{T{0}}, den{T{0}};
|
|
|
|
for (int64_t i = 0; i < n; ++i){
|
|
T h = x - xx[i];
|
|
if (h == T{0}){
|
|
return yy[i];
|
|
}else{
|
|
T temp = w[i]/h;
|
|
num += temp*yy[i];
|
|
den += temp;
|
|
}
|
|
}
|
|
return num/den;
|
|
}
|
|
|
|
|
|
T interp(T x) {
|
|
return rawinterp(1, x);
|
|
}
|
|
|
|
};
|
|
|
|
} // namespace numerics
|
|
|