Skip to content

Commit

Permalink
added some test
Browse files Browse the repository at this point in the history
  • Loading branch information
@BryanFlynt
BryanFlynt committed Sep 3, 2021
1 parent f33fcdc commit 2605554
Show file tree
Hide file tree
Showing 8 changed files with 389 additions and 88 deletions.
91 changes: 49 additions & 42 deletions include/polycalc/interpolation/lagrange.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#include <limits>
#include <vector>

#include "parameters.hpp"
#include "polycalc/parameters.hpp"

namespace polycalc {
namespace interpolation {
Expand All @@ -37,20 +37,16 @@ class Lagrange {
*
* @param points Iterable container of point locations
*/
template <IterableContainer>
template <typename IterableContainer>
Lagrange(const IterableContainer& points);

/**
* Assign the nodal locations of interpolation
*
* @param points Iterable container of point locations
*/
template <IterableContainer>
void set_points(const IterableContainer& points) {
x_.clear();
std::copy(points.begin(), points.end(), std::back_inserter(x_));
this->init_();
}
template <typename IterableContainer>
void reset(const IterableContainer& points);

/**
* Return the i'th polynomial evaluated at location x
Expand Down Expand Up @@ -97,7 +93,7 @@ class Lagrange {

protected:
/// Precompute the denominators
void init_();
void precompute_();

private:
std::vector<value_type> x_; ///< Node Locations
Expand All @@ -108,55 +104,67 @@ namespace detail {

/// Perform Subtraction with less Subtractive Cancellation
/**
*
* When values are close it strips off the simular parts and scales up the
* remainder to perform the subtraction. Then unscales the difference of the
* remainder to proces a "better" answer.
*
* @note
* Expensive so only use when necessary.
*/
template <typename T>
T calc_diff(const T x, const T y) {
// template <typename T>
// T calc_diff(const T x, const T y) {

// If the values are "far" away then skip
if (std::abs(x - y) < static_cast<T>(0.0001)) {
return (x - y);
}
// // If the values are "far" away then skip
// if ( static_cast<T>(1.0E-4) < std::abs(x - y) ){
// return (x - y);
// }

using int_type = long int;
using real_type = long double;
// using int_type = long int;
// using real_type = long double;

constexpr auto eps = std::numeric_limits<real_type>::epsilon();
const auto digits = static_cast<int_type>(std::abs(std::floor(std::log10(eps))) - 1);
const auto sub_expn = static_cast<int_type>(std::floor(std::log10(std::abs(x - y))));
assert(sub_expn < 0);
// constexpr auto eps = std::numeric_limits<real_type>::epsilon();
// const auto digits = static_cast<int_type>(std::abs(std::floor(std::log10(eps))) - 1);
// const auto sub_expn = static_cast<int_type>(std::floor(std::log10(std::abs(x - y))));
// assert(sub_expn < 0);

const real_type scaled_x = x * std::pow(10ULL, -sub_expn); // Scale up by non-overlaping digits
const real_type remain_x = scaled_x - std::floor(scaled_x); // Get remainder that is different
const real_type scaled_rx = remain_x * std::pow(10ULL, digits+sub_expn); // Scale up remainder to max precision
const real_type trimed_rx = std::round(scaled_rx); // Remove junk that cannot be represented
// const real_type scaled_x = x * std::pow(10ULL, -sub_expn); // Scale up by non-overlaping digits
// const real_type remain_x = scaled_x - std::floor(scaled_x); // Get remainder that is different
// const real_type scaled_rx = remain_x * std::pow(10ULL, digits+sub_expn); // Scale up remainder to max precision
// const real_type trimed_rx = std::round(scaled_rx); // Remove junk that cannot be represented

const real_type scaled_y = y * std::pow(10ULL, -sub_expn); // Scale up by non-overlaping digits
const real_type remain_y = scaled_y - std::floor(scaled_y); // Get remainder that is different
const real_type scaled_ry = remain_y * std::pow(10ULL, digits+sub_expn); // Scale up remainder to may precision
const real_type trimed_ry = std::round(scaled_ry); // Remove junk that cannot be represented
// const real_type scaled_y = y * std::pow(10ULL, -sub_expn); // Scale up by non-overlaping digits
// const real_type remain_y = scaled_y - std::floor(scaled_y); // Get remainder that is different
// const real_type scaled_ry = remain_y * std::pow(10ULL, digits+sub_expn); // Scale up remainder to may precision
// const real_type trimed_ry = std::round(scaled_ry); // Remove junk that cannot be represented

return static_cast<T>(trimed_rx - trimed_ry) * std::pow(10ULL, -digits); // Unscale the difference of remainder
// return static_cast<T>(trimed_rx - trimed_ry) * std::pow(10ULL, -digits); // Unscale the difference of remainder
// };


template <typename T>
T calc_diff(const T x, const T y) {
return (x - y);
};

} /* namespace detail */

template <typename T, typename P>
template <IterableContainer>
template <typename IterableContainer>
Lagrange<T, P>::Lagrange(const IterableContainer& points) {
this->set_points(points);
this->reset(points);
}

template <typename T, typename P>
template <IterableContainer>
void Lagrange<T, P>::set_points(const IterableContainer& points) {
template <typename IterableContainer>
void Lagrange<T, P>::reset(const IterableContainer& points) {
x_.clear();
std::copy(points.begin(), points.end(), std::back_inserter(x_));
this->init_();
this->precompute_();
}

template <typename T, typename P>
typename Lagrange<T, P>::value_type Lagrange<T, P>::eval(const size_type i, const value_type x) const {
assert(i >= 0);
assert(i < x_.size());

value_type product = 1;
Expand All @@ -172,8 +180,8 @@ typename Lagrange<T, P>::value_type Lagrange<T, P>::eval(const size_type i, cons

template <typename T, typename P>
typename Lagrange<T, P>::value_type Lagrange<T, P>::ddx(const size_type i, const value_type x) const {
assert(i >= 0);
assert(i < x_.size());
using std::abs;

const value_type one = 1;

Expand All @@ -184,7 +192,7 @@ typename Lagrange<T, P>::value_type Lagrange<T, P>::ddx(const size_type i, const
for (size_type j = 0; j < x_.size(); ++j) {
if (i != j) {
dist = detail::calc_diff(x, x_[j]);
if (abs(dist) < params::TOL) {
if (std::abs(dist) < params::TOL) {
x_equals_point = true;
} else {
rsum += one / dist;
Expand All @@ -200,12 +208,11 @@ typename Lagrange<T, P>::value_type Lagrange<T, P>::ddx(const size_type i, const
}

template <typename T, typename P>
void Lagrange<T, P>::init_() {
using std::abs;
void Lagrange<T, P>::precompute_() {
const size_type sz = x_.size();

d_.resize(sz);
std::fill(d_.begin().d_.end(), 1);
std::fill(d_.begin(),d_.end(), 1);
for (size_type j = 0; j < sz; ++j) {
for (size_type i = 0; i < j; ++i) {
d_[j] *= detail::calc_diff(x_[j], x_[i]);
Expand All @@ -214,7 +221,7 @@ void Lagrange<T, P>::init_() {
d_[j] *= detail::calc_diff(x_[j], x_[i]);
}
// Test for distinct points
assert(abs(d_[j]) > params::TOL);
assert(std::abs(d_[j]) > params::TOL);
}
}

Expand Down
6 changes: 2 additions & 4 deletions include/polycalc/quadrature/gauss_jacobi.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,6 @@ std::vector<typename GaussJacobi<T, P>::value_type> GaussJacobi<T, P>::zeros(con
template <typename T, typename P>
std::vector<typename GaussJacobi<T, P>::value_type> GaussJacobi<T, P>::weights(const unsigned n) const {
assert(n > 0);
using std::pow;
using std::tgamma;

// Get location of zeros
auto z = this->zeros(n);
Expand All @@ -75,8 +73,8 @@ std::vector<typename GaussJacobi<T, P>::value_type> GaussJacobi<T, P>::weights(c
const value_type np1 = 1 + n;

value_type fac;
fac = pow(two, apb + one) * tgamma(poly_.alpha + np1) * tgamma(poly_.beta + np1);
fac /= tgamma(np1) * tgamma(apb + np1);
fac = std::pow(two, apb + one) * std::tgamma(poly_.alpha + np1) * std::tgamma(poly_.beta + np1);
fac /= std::tgamma(np1) * std::tgamma(apb + np1);

for (size_type i = 0; i < n; ++i) {
w[i] = fac / (w[i] * w[i] * (one - z[i] * z[i]));
Expand Down
6 changes: 2 additions & 4 deletions include/polycalc/quadrature/gauss_lobatto.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,6 @@ std::vector<typename GaussLobatto<T, P>::value_type> GaussLobatto<T, P>::zeros(c
template <typename T, typename P>
std::vector<typename GaussLobatto<T, P>::value_type> GaussLobatto<T, P>::weights(const unsigned n) const {
assert(n > 0);
using std::pow;
using std::tgamma;

// Weights to return
std::vector<value_type> w(n);
Expand All @@ -102,8 +100,8 @@ std::vector<typename GaussLobatto<T, P>::value_type> GaussLobatto<T, P>::weights
const value_type apb = alpha_ + beta_;

value_type fac;
fac = pow(two, apb + one) * tgamma(alpha_ + n) * tgamma(beta_ + n);
fac /= (n - 1) * tgamma(n) * tgamma(alpha_ + beta_ + n + one);
fac = std::pow(two, apb + one) * std::tgamma(alpha_ + n) * std::tgamma(beta_ + n);
fac /= (n - 1) * std::tgamma(n) * std::tgamma(alpha_ + beta_ + n + one);

for (size_type i = 0; i < n; ++i) {
w[i] = fac / (w[i] * w[i]);
Expand Down
10 changes: 4 additions & 6 deletions include/polycalc/quadrature/gauss_radau.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,6 @@ template <typename T, typename P>
template <int Z1>
std::vector<typename GaussRadau<T, P>::value_type> GaussRadau<T, P>::weights(const unsigned n) const {
assert(n > 0);
using std::pow;
using std::tgamma;

if (n == 1) {
return std::vector<value_type>(1, 2);
Expand All @@ -106,16 +104,16 @@ std::vector<typename GaussRadau<T, P>::value_type> GaussRadau<T, P>::weights(con
value_type fac;

if constexpr (Z1 == +1) { // Z = +1
fac = pow(two, apb) * tgamma(alpha_ + n) * tgamma(beta_ + n);
fac /= tgamma(n) * (alpha_ + n) * tgamma(apb + n + 1);
fac = std::pow(two, apb) * std::tgamma(alpha_ + n) * std::tgamma(beta_ + n);
fac /= std::tgamma(n) * (alpha_ + n) * std::tgamma(apb + n + 1);

for (size_type i = 0; i < n; ++i) {
w[i] = fac * (one + z[i]) / (w[i] * w[i]);
}
w[n - 1] *= (alpha_ + one);
} else if constexpr (Z1 == -1) { // Z = -1
fac = pow(two, apb) * tgamma(alpha_ + n) * tgamma(beta_ + n);
fac /= tgamma(n) * (beta_ + n) * tgamma(apb + n + 1);
fac = std::pow(two, apb) * std::tgamma(alpha_ + n) * std::tgamma(beta_ + n);
fac /= std::tgamma(n) * (beta_ + n) * std::tgamma(apb + n + 1);

for (size_type i = 0; i < n; ++i) {
w[i] = fac * (one - z[i]) / (w[i] * w[i]);
Expand Down
3 changes: 2 additions & 1 deletion test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,4 +43,5 @@ function(add_catch_test test_name)
endfunction()


add_catch_test(test_jacobi)
add_catch_test(jacobi)
add_catch_test(lagrange)
Loading

0 comments on commit 2605554

Please sign in to comment.