From cfd5a353f63f977ee13ef4bbdccea36bcdf36348 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Mon, 15 Jul 2019 10:42:49 -0400 Subject: [PATCH] Cardinal B-splines --- doc/math.qbk | 2 + test/Jamfile.v2 | 1 + test/cardinal_b_spline_test.cpp | 197 ++++++++++++++++++++++++++++++++ 3 files changed, 200 insertions(+) create mode 100644 test/cardinal_b_spline_test.cpp diff --git a/doc/math.qbk b/doc/math.qbk index 3b2723dfd..8373c5780 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -300,6 +300,7 @@ and use the function's name as the link text.] [def __legendre_q [link math_toolkit.sf_poly.legendre legendre_q]] [def __legendre_next [link math_toolkit.sf_poly.legendre legendre_next]] [def __hermite [link math_toolkit.sf_poly.hermite hermite]] +[def __cardinal_b_splines [link math_toolkit.sf_poly.cardinal_b_splines cardinal_b_splines]] [/Misc] [def __expint [link math_toolkit.expint.expint_i expint]] @@ -600,6 +601,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include sf/hermite.qbk] [include sf/chebyshev.qbk] [include sf/spherical_harmonic.qbk] +[include sf/cardinal_b_splines.qbk] [endsect] [/section:sf_poly Polynomials] [section:bessel Bessel Functions] diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 53c7a9e64..5203255ab 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -902,6 +902,7 @@ test-suite misc : [ run test_vector_barycentric_rational.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] [ run test_constant_generate.cpp : : : release USE_CPP_FLOAT=1 off:no ] [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] off msvc:/bigobj release ] + [ run cardinal_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run whittaker_shannon_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] ] [ run cardinal_quadratic_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ] diff --git a/test/cardinal_b_spline_test.cpp b/test/cardinal_b_spline_test.cpp new file mode 100644 index 000000000..89a63827a --- /dev/null +++ b/test/cardinal_b_spline_test.cpp @@ -0,0 +1,197 @@ +/* + * Copyright Nick Thompson, 2019 + * 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 +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + +using std::abs; +using boost::math::cardinal_b_spline; +using boost::math::forward_cardinal_b_spline; + +template +void test_box() +{ + Real t = cardinal_b_spline<0>(1.1); + Real expected = 0; + CHECK_ULP_CLOSE(expected, t, 0); + + t = cardinal_b_spline<0>(-1.1); + expected = 0; + CHECK_ULP_CLOSE(expected, t, 0); + + Real h = Real(1)/Real(256); + for (Real t = -Real(1)/Real(2)+h; t < Real(1)/Real(2); t += h) + { + expected = 1; + CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0); + } + + for (Real t = h; t < 1; t += h) + { + expected = 1; + CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<0>(t), 0); + } +} + +template +void test_hat() +{ + Real t = cardinal_b_spline<1>(2.1); + Real expected = 0; + CHECK_ULP_CLOSE(expected, t, 0); + + t = cardinal_b_spline<1>(-2.1); + expected = 0; + CHECK_ULP_CLOSE(expected, t, 0); + + Real h = Real(1)/Real(256); + for (Real t = -1; t <= 1; t += h) + { + expected = 1-abs(t); + if(!CHECK_ULP_CLOSE(expected, cardinal_b_spline<1>(t), 0) ) + { + std::cerr << " Problem at t = " << t << "\n"; + } + } + + for (Real t = 0; t < 2; t += h) + { + expected = 1 - abs(t-1); + CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<1>(t), 0); + } +} + +template +void test_quadratic() +{ + using std::abs; + auto b2 = [](Real x) { + Real absx = abs(x); + if (absx >= 3/Real(2)) { + return Real(0); + } + if (absx >= 1/Real(2)) { + Real t = absx - 3/Real(2); + return t*t/2; + } + Real t1 = absx - 1/Real(2); + Real t2 = absx + 1/Real(2); + return (2-t1*t1 -t2*t2)/2; + }; + + Real h = 1/Real(256); + for (Real t = -5; t <= 5; t += h) { + Real expected = b2(t); + CHECK_ULP_CLOSE(expected, cardinal_b_spline<2>(t), 0); + } +} + +template +void test_cubic() +{ + Real expected = Real(2)/Real(3); + Real computed = cardinal_b_spline<3, Real>(0); + CHECK_ULP_CLOSE(expected, computed, 0); + + expected = Real(1)/Real(6); + computed = cardinal_b_spline<3, Real>(1); + CHECK_ULP_CLOSE(expected, computed, 0); + + expected = Real(0); + computed = cardinal_b_spline<3, Real>(2); + CHECK_ULP_CLOSE(expected, computed, 0); +} + +template +void test_quintic() +{ + Real expected = Real(11)/Real(20); + Real computed = cardinal_b_spline<5, Real>(0); + CHECK_ULP_CLOSE(expected, computed, 0); + + expected = Real(13)/Real(60); + computed = cardinal_b_spline<5, Real>(1); + CHECK_ULP_CLOSE(expected, computed, 1); + + expected = Real(1)/Real(120); + computed = cardinal_b_spline<5, Real>(2); + CHECK_ULP_CLOSE(expected, computed, 0); + + expected = Real(0); + computed = cardinal_b_spline<5, Real>(3); + CHECK_ULP_CLOSE(expected, computed, 0); + +} + +template +void test_partition_of_unity() +{ + std::mt19937 gen(323723); + Real supp = (n+1.0)/2.0; + std::uniform_real_distribution dis(-supp, -supp+1); + + for(size_t i = 0; i < 500; ++i) { + Real x = dis(gen); + Real one = 0; + while (x < supp) { + one += cardinal_b_spline(x); + x += 1; + } + if(!CHECK_ULP_CLOSE(1, one, n)) { + std::cerr << " Partition of unity failure at n = " << n << "\n"; + } + } +} + + +int main() +{ + test_box(); + test_box(); + test_box(); + + test_hat(); + test_hat(); + test_hat(); + + test_quadratic(); + test_quadratic(); + test_quadratic(); + + test_cubic(); + test_cubic(); + test_cubic(); + + test_quintic(); + test_quintic(); + test_quintic(); + + test_partition_of_unity<1, double>(); + test_partition_of_unity<2, double>(); + test_partition_of_unity<3, double>(); + test_partition_of_unity<4, double>(); + test_partition_of_unity<5, double>(); + test_partition_of_unity<6, double>(); + +#ifdef BOOST_HAS_FLOAT128 + test_box(); + test_hat(); + test_quadratic(); + test_cubic(); + test_quintic(); +#endif + + return boost::math::test::report_errors(); +}