generic-library/boost
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
boost/libs/math/test/cubic_hermite_test.cpp
Edouard DUPIN d0115b733d [DEV] add v1.76.0
2021-10-05 21:37:46 +02:00

462 lines
14 KiB
C++
Raw Permalink Blame History

/*
* Copyright Nick Thompson, 2020
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include "math_unit_test.hpp"
#include <numeric>
#include <utility>
#include <random>
#include <array>
#include <vector>
#include <boost/math/interpolators/cubic_hermite.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/circular_buffer.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif
using boost::math::interpolators::cubic_hermite;
using boost::math::interpolators::cardinal_cubic_hermite;
using boost::math::interpolators::cardinal_cubic_hermite_aos;
template<typename Real>
void test_constant()
{
Real x0 = 0;
std::vector<Real> x{x0,1,2,3, 9, 22, 81};
std::vector<Real> y(x.size());
for (auto & t : y)
{
t = 7;
}
std::vector<Real> dydx(x.size(), Real(0));
auto x_copy = x;
auto y_copy = y;
auto dydx_copy = dydx;
auto hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
// Now check the boundaries:
Real tlo = x.front();
Real thi = x.back();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2);
CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
boost::circular_buffer<Real> x_buf(x.size());
for (auto & t : x) {
x_buf.push_back(t);
}
boost::circular_buffer<Real> y_buf(x.size());
for (auto & t : y) {
y_buf.push_back(t);
}
boost::circular_buffer<Real> dydx_buf(x.size());
for (auto & t : dydx) {
dydx_buf.push_back(t);
}
auto circular_hermite_spline = cubic_hermite(std::move(x_buf), std::move(y_buf), std::move(dydx_buf));
for (Real t = x[0]; t <= x.back(); t += 0.25) {
CHECK_ULP_CLOSE(Real(7), circular_hermite_spline(t), 2);
CHECK_ULP_CLOSE(Real(0), circular_hermite_spline.prime(t), 2);
}
circular_hermite_spline.push_back(x.back() + 1, 7, 0);
CHECK_ULP_CLOSE(Real(0), circular_hermite_spline.prime(x.back()+1), 2);
}
template<typename Real>
void test_linear()
{
std::vector<Real> x{0,1,2,3};
std::vector<Real> y{0,1,2,3};
std::vector<Real> dydx{1,1,1,1};
auto x_copy = x;
auto y_copy = y;
auto dydx_copy = dydx;
auto hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
CHECK_ULP_CLOSE(y[0], hermite_spline(x[0]), 0);
CHECK_ULP_CLOSE(Real(1)/Real(2), hermite_spline(Real(1)/Real(2)), 10);
CHECK_ULP_CLOSE(y[1], hermite_spline(x[1]), 0);
CHECK_ULP_CLOSE(Real(3)/Real(2), hermite_spline(Real(3)/Real(2)), 10);
CHECK_ULP_CLOSE(y[2], hermite_spline(x[2]), 0);
CHECK_ULP_CLOSE(Real(5)/Real(2), hermite_spline(Real(5)/Real(2)), 10);
CHECK_ULP_CLOSE(y[3], hermite_spline(x[3]), 0);
x.resize(45);
y.resize(45);
dydx.resize(45);
for (size_t i = 0; i < x.size(); ++i) {
x[i] = i;
y[i] = i;
dydx[i] = 1;
}
x_copy = x;
y_copy = y;
dydx_copy = dydx;
hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
for (Real t = 0; t < x.back(); t += 0.5) {
CHECK_ULP_CLOSE(t, hermite_spline(t), 0);
CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(t), 0);
}
boost::circular_buffer<Real> x_buf(x.size());
for (auto & t : x) {
x_buf.push_back(t);
}
boost::circular_buffer<Real> y_buf(x.size());
for (auto & t : y) {
y_buf.push_back(t);
}
boost::circular_buffer<Real> dydx_buf(x.size());
for (auto & t : dydx) {
dydx_buf.push_back(t);
}
auto circular_hermite_spline = cubic_hermite(std::move(x_buf), std::move(y_buf), std::move(dydx_buf));
for (Real t = x[0]; t <= x.back(); t += 0.25) {
CHECK_ULP_CLOSE(t, circular_hermite_spline(t), 2);
CHECK_ULP_CLOSE(Real(1), circular_hermite_spline.prime(t), 2);
}
circular_hermite_spline.push_back(x.back() + 1, y.back()+1, 1);
CHECK_ULP_CLOSE(Real(y.back() + 1), circular_hermite_spline(Real(x.back()+1)), 2);
CHECK_ULP_CLOSE(Real(1), circular_hermite_spline.prime(Real(x.back()+1)), 2);
}
template<typename Real>
void test_quadratic()
{
std::vector<Real> x(50);
std::default_random_engine rd;
std::uniform_real_distribution<Real> dis(0.1,1);
Real x0 = dis(rd);
x[0] = x0;
for (size_t i = 1; i < x.size(); ++i) {
x[i] = x[i-1] + dis(rd);
}
Real xmax = x.back();
std::vector<Real> y(x.size());
std::vector<Real> dydx(x.size());
for (size_t i = 0; i < x.size(); ++i) {
y[i] = x[i]*x[i]/2;
dydx[i] = x[i];
}
auto s = cubic_hermite(std::move(x), std::move(y), std::move(dydx));
for (Real t = x0; t <= xmax; t+= 0.0125)
{
CHECK_ULP_CLOSE(t*t/2, s(t), 5);
CHECK_ULP_CLOSE(t, s.prime(t), 138);
}
}
template<typename Real>
void test_interpolation_condition()
{
for (size_t n = 4; n < 50; ++n) {
std::vector<Real> x(n);
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::default_random_engine rd;
std::uniform_real_distribution<Real> dis(0,1);
Real x0 = dis(rd);
x[0] = x0;
y[0] = dis(rd);
for (size_t i = 1; i < n; ++i) {
x[i] = x[i-1] + dis(rd);
y[i] = dis(rd);
dydx[i] = dis(rd);
}
auto x_copy = x;
auto y_copy = y;
auto dydx_copy = dydx;
auto s = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
//std::cout << "s = " << s << "\n";
for (size_t i = 0; i < x.size(); ++i) {
CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2);
}
}
}
template<typename Real>
void test_cardinal_constant()
{
Real x0 = 0;
Real dx = 2;
std::vector<Real> y(25);
for (auto & t : y) {
t = 7;
}
std::vector<Real> dydx(y.size(), Real(0));
auto hermite_spline = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
for (Real t = x0; t <= x0 + 24*dx; t += 0.25) {
CHECK_ULP_CLOSE(Real(7), hermite_spline(t), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(t), 2);
}
// Array of structs:
std::vector<std::array<Real, 2>> data(25);
for (auto & t : data) {
t[0] = 7;
t[1] = 0;
}
auto hermite_spline_aos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
for (Real t = x0; t <= x0 + 24*dx; t += 0.25) {
if (!CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(t), 2)) {
std::cerr << " Wrong evaluation at t = " << t << "\n";
}
if (!CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(t), 2)) {
std::cerr << " Wrong evaluation at t = " << t << "\n";
}
}
// Now check the boundaries:
Real tlo = x0;
Real thi = x0 + (25-1)*dx;
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2);
CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2);
CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(thi), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(thi), 2);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>
void test_cardinal_linear()
{
Real x0 = 0;
Real dx = 1;
std::vector<Real> y{0,1,2,3};
std::vector<Real> dydx{1,1,1,1};
auto y_copy = y;
auto dydx_copy = dydx;
auto hermite_spline = cardinal_cubic_hermite(std::move(y_copy), std::move(dydx_copy), x0, dx);
CHECK_ULP_CLOSE(y[0], hermite_spline(0), 0);
CHECK_ULP_CLOSE(Real(1)/Real(2), hermite_spline(Real(1)/Real(2)), 10);
CHECK_ULP_CLOSE(y[1], hermite_spline(1), 0);
CHECK_ULP_CLOSE(Real(3)/Real(2), hermite_spline(Real(3)/Real(2)), 10);
CHECK_ULP_CLOSE(y[2], hermite_spline(2), 0);
CHECK_ULP_CLOSE(Real(5)/Real(2), hermite_spline(Real(5)/Real(2)), 10);
CHECK_ULP_CLOSE(y[3], hermite_spline(3), 0);
y.resize(45);
dydx.resize(45);
for (size_t i = 0; i < y.size(); ++i) {
y[i] = i;
dydx[i] = 1;
}
hermite_spline = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
for (Real t = 0; t < 44; t += 0.5) {
CHECK_ULP_CLOSE(t, hermite_spline(t), 0);
CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(t), 0);
}
std::vector<std::array<Real, 2>> data(45);
for (size_t i = 0; i < data.size(); ++i) {
data[i][0] = i;
data[i][1] = 1;
}
auto hermite_spline_aos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
for (Real t = 0; t < 44; t += 0.5) {
CHECK_ULP_CLOSE(t, hermite_spline_aos(t), 0);
CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(t), 0);
}
Real tlo = x0;
Real thi = x0 + (45-1)*dx;
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(tlo), hermite_spline(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), hermite_spline(thi), 2);
CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(thi), 2);
CHECK_ULP_CLOSE(Real(tlo), hermite_spline_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), hermite_spline_aos(thi), 2);
CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(thi), 2);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>
void test_cardinal_quadratic()
{
Real x0 = -1;
Real dx = Real(1)/Real(256);
std::vector<Real> y(50);
std::vector<Real> dydx(y.size());
for (size_t i = 0; i < y.size(); ++i) {
Real x = x0 + i*dx;
y[i] = x*x/2;
dydx[i] = x;
}
auto s = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
for (Real t = x0; t <= x0 + 49*dx; t+= 0.0125)
{
CHECK_ULP_CLOSE(t*t/2, s(t), 12);
CHECK_ULP_CLOSE(t, s.prime(t), 70);
}
std::vector<std::array<Real, 2>> data(50);
for (size_t i = 0; i < data.size(); ++i) {
Real x = x0 + i*dx;
data[i][0] = x*x/2;
data[i][1] = x;
}
auto saos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
for (Real t = x0; t <= x0 + 49*dx; t+= 0.0125)
{
CHECK_ULP_CLOSE(t*t/2, saos(t), 12);
CHECK_ULP_CLOSE(t, saos.prime(t), 70);
}
auto [tlo, thi] = s.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(tlo*tlo/2), s(tlo), 3);
CHECK_ULP_CLOSE(Real(thi*thi/2), s(thi), 3);
CHECK_ULP_CLOSE(Real(tlo), s.prime(tlo), 3);
CHECK_ULP_CLOSE(Real(thi), s.prime(thi), 3);
CHECK_ULP_CLOSE(Real(tlo*tlo/2), saos(tlo), 3);
CHECK_ULP_CLOSE(Real(thi*thi/2), saos(thi), 3);
CHECK_ULP_CLOSE(Real(tlo), saos.prime(tlo), 3);
CHECK_ULP_CLOSE(Real(thi), saos.prime(thi), 3);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>
void test_cardinal_interpolation_condition()
{
for (size_t n = 4; n < 50; ++n) {
std::vector<Real> y(n);
std::vector<Real> dydx(n);
std::default_random_engine rd;
std::uniform_real_distribution<Real> dis(0.1,1);
Real x0 = Real(2);
Real dx = Real(1)/Real(128);
for (size_t i = 0; i < n; ++i) {
y[i] = dis(rd);
dydx[i] = dis(rd);
}
auto y_copy = y;
auto dydx_copy = dydx;
auto s = cardinal_cubic_hermite(std::move(y_copy), std::move(dydx_copy), x0, dx);
for (size_t i = 0; i < y.size(); ++i) {
CHECK_ULP_CLOSE(y[i], s(x0 + i*dx), 2);
CHECK_ULP_CLOSE(dydx[i], s.prime(x0 + i*dx), 2);
}
}
}
int main()
{
test_constant<float>();
test_linear<float>();
test_quadratic<float>();
test_interpolation_condition<float>();
test_cardinal_constant<float>();
test_cardinal_linear<float>();
test_cardinal_quadratic<float>();
test_cardinal_interpolation_condition<float>();
test_constant<double>();
test_linear<double>();
test_quadratic<double>();
test_interpolation_condition<double>();
test_cardinal_constant<double>();
test_cardinal_linear<double>();
test_cardinal_quadratic<double>();
test_cardinal_interpolation_condition<double>();
test_constant<long double>();
test_linear<long double>();
test_quadratic<long double>();
test_interpolation_condition<long double>();
test_cardinal_constant<long double>();
test_cardinal_linear<long double>();
test_cardinal_quadratic<long double>();
test_cardinal_interpolation_condition<long double>();
#ifdef BOOST_HAS_FLOAT128
test_constant<float128>();
test_linear<float128>();
test_cardinal_constant<float128>();
test_cardinal_linear<float128>();
#endif
return boost::math::test::report_errors();
}
Reference in New Issue View Git Blame Copy Permalink