From e327be5887a5fd8aa4ae3544b6738bdb50d43f10 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 11 Aug 2019 14:52:47 -0400 Subject: [PATCH 1/4] Cardinal B-spline derivatives [CI SKIP] --- doc/sf/cardinal_b_splines.qbk | 14 +- .../special_functions/cardinal_b_spline.hpp | 172 +++++++++++++----- test/cardinal_b_spline_test.cpp | 13 ++ 3 files changed, 148 insertions(+), 51 deletions(-) diff --git a/doc/sf/cardinal_b_splines.qbk b/doc/sf/cardinal_b_splines.qbk index 7f7c4e15e..f4e2f3938 100644 --- a/doc/sf/cardinal_b_splines.qbk +++ b/doc/sf/cardinal_b_splines.qbk @@ -16,7 +16,13 @@ namespace boost{ namespace math{ template - Real cardinal_b_spline(Real x); + auto cardinal_b_spline(Real x); + + template + auto cardinal_b_spline_prime(Real x); + + template + auto cardinal_b_spline_double_prime(Real x); template Real forward_cardinal_b_spline(Real x); @@ -34,11 +40,17 @@ For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere A basic usage is as follows: using boost::math::cardinal_b_spline; + using boost::math::cardinal_b_spline_prime; + using boost::math::cardinal_b_spline_double_prime; double x = 0.5; // B₀(x), the box function: double y = cardinal_b_spline<0>(x); // B₁(x), the hat function: y = cardinal_b_spline<1>(x); + // First derivative of B₃: + yp = cardinal_b_spline_prime<3>(x); + // Second derivative of B₃: + ypp = cardinal_b_spline_double_prime<3>(x); [$../graphs/central_b_splines.svg] diff --git a/include/boost/math/special_functions/cardinal_b_spline.hpp b/include/boost/math/special_functions/cardinal_b_spline.hpp index 399934ffd..4e40058ed 100644 --- a/include/boost/math/special_functions/cardinal_b_spline.hpp +++ b/include/boost/math/special_functions/cardinal_b_spline.hpp @@ -8,6 +8,8 @@ #include #include +#include +#include namespace boost { namespace math { @@ -30,65 +32,135 @@ namespace detail { template Real cardinal_b_spline(Real x) { - if (x < 0) { - // All B-splines are even functions: - return cardinal_b_spline(-x); - } + static_assert(!std::is_integral::value, "Does not work with integral types."); - // Someday we'll be able to 'if constexpr' on many compilers. - // Not today though! In any case, it's not a big performance hit. - if /*constexpr*/ (n==0) - { - if (x < Real(1)/Real(2)) { - return Real(1); + if (x < 0) { + // All B-splines are even functions: + return cardinal_b_spline(-x); } - else if (x == Real(1)/Real(2)) { - return Real(1)/Real(2); - } - else { - return Real(0); - } - } - if /*constexpr*/ (n==1) { - return detail::B1(x); - } - - Real supp_max = (n+1)/Real(2); - if (x >= supp_max) - { - return Real(0); - } - - // Fill v with values of B1: - // At most two of these terms are nonzero, and at least 1. - // There is only one non-zero term when n is odd and x = 0. - std::array v; - Real z = x + 1 - supp_max; - for (unsigned i = 0; i < n; ++i) - { - v[i] = detail::B1(z); - z += 1; - } - - Real smx = supp_max - x; - for (unsigned j = 2; j <= n; ++j) - { - Real a = (j + 1 - smx); - Real b = smx; - for(unsigned k = 0; k <= n - j; ++k) + if (n==0) { - v[k] = (a*v[k+1] + b*v[k])/Real(j); - a += 1; - b -= 1; + if (x < Real(1)/Real(2)) { + return Real(1); + } + else if (x == Real(1)/Real(2)) { + return Real(1)/Real(2); + } + else { + return Real(0); + } } - } - return v[0]; + if (n==1) + { + return detail::B1(x); + } + + Real supp_max = (n+1)/Real(2); + if (x >= supp_max) + { + return Real(0); + } + + // Fill v with values of B1: + // At most two of these terms are nonzero, and at least 1. + // There is only one non-zero term when n is odd and x = 0. + std::array v; + Real z = x + 1 - supp_max; + for (unsigned i = 0; i < n; ++i) + { + v[i] = detail::B1(z); + z += 1; + } + + Real smx = supp_max - x; + for (unsigned j = 2; j <= n; ++j) + { + Real a = (j + 1 - smx); + Real b = smx; + for(unsigned k = 0; k <= n - j; ++k) + { + v[k] = (a*v[k+1] + b*v[k])/Real(j); + a += 1; + b -= 1; + } + } + + return v[0]; +} + + +template +Real cardinal_b_spline_prime(Real x) +{ + static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integer types."); + + if (x < 0) + { + // All B-splines are even functions, so derivatives are odd: + return -cardinal_b_spline_prime(-x); + } + + + if (n==0) + { + // Kinda crazy but you get what you ask for! + if (x == Real(1)/Real(2)) + { + return std::numeric_limits::infinity(); + } + else + { + return Real(0); + } + } + + Real supp_max = (n+1)/Real(2); + if (x >= supp_max) + { + return Real(0); + } + + if (n==1) { + if (x==0) { + return Real(0); + } + return Real(-1); + } + + // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): + std::array v; + Real z = x + 1 - supp_max; + for (unsigned i = 0; i < n; ++i) + { + v[i] = detail::B1(z); + z += 1; + } + + Real smx = supp_max - x; + for (unsigned j = 2; j <= n - 1; ++j) + { + Real a = (j + 1 - smx); + Real b = smx; + for(unsigned k = 0; k <= n - j; ++k) + { + v[k] = (a*v[k+1] + b*v[k])/Real(j); + a += 1; + b -= 1; + } + } + + return v[1] - v[0]; } template -Real forward_cardinal_b_spline(Real x) { +auto forward_cardinal_b_spline(Real x) +{ + if constexpr (std::is_integral::value) + { + return forward_cardinal_b_spline(x); + } return cardinal_b_spline(x - (n+1)/Real(2)); } diff --git a/test/cardinal_b_spline_test.cpp b/test/cardinal_b_spline_test.cpp index b819abba6..20c5c5dc7 100644 --- a/test/cardinal_b_spline_test.cpp +++ b/test/cardinal_b_spline_test.cpp @@ -19,6 +19,7 @@ using boost::multiprecision::float128; using std::abs; using boost::math::cardinal_b_spline; +using boost::math::cardinal_b_spline_prime; using boost::math::forward_cardinal_b_spline; template @@ -27,6 +28,7 @@ void test_box() Real t = cardinal_b_spline<0>(Real(1.1)); Real expected = 0; CHECK_ULP_CLOSE(expected, t, 0); + CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0); t = cardinal_b_spline<0>(Real(-1.1)); expected = 0; @@ -37,6 +39,8 @@ void test_box() { expected = 1; CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0); + expected = 0; + CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0); } for (t = h; t < 1; t += h) @@ -65,6 +69,15 @@ void test_hat() { std::cerr << " Problem at t = " << t << "\n"; } + if (t < 0) { + CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0); + } + else if (t == 0) { + CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime<1>(t), 0); + } + else if (t > 0) { + CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime<1>(t), 0); + } } for (t = 0; t < 2; t += h) From 5c8fcb4cd592707f3910093e3c4694c371e1f1e5 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Mon, 12 Aug 2019 10:07:06 -0400 Subject: [PATCH 2/4] Cardinal B-spline derivatives. --- doc/graphs/central_b_spline_derivatives.svg | 56 +++++++++++++ .../central_b_spline_second_derivatives.svg | 54 ++++++++++++ doc/sf/cardinal_b_splines.qbk | 2 + .../detail/cubic_b_spline_detail.hpp | 2 +- .../special_functions/cardinal_b_spline.hpp | 72 +++++++++++++--- test/cardinal_b_spline_test.cpp | 83 ++++++++++++++++++- 6 files changed, 255 insertions(+), 14 deletions(-) create mode 100644 doc/graphs/central_b_spline_derivatives.svg create mode 100644 doc/graphs/central_b_spline_second_derivatives.svg diff --git a/doc/graphs/central_b_spline_derivatives.svg b/doc/graphs/central_b_spline_derivatives.svg new file mode 100644 index 000000000..490f93fd5 --- /dev/null +++ b/doc/graphs/central_b_spline_derivatives.svg @@ -0,0 +1,56 @@ + + + +B-spline Derivatives + + + + +-0.75 + +-0.5 + +-0.25 + +0 + +0.25 + +0.5 + +0.75 + +1 + +-1.6 + +-1.2 + +-0.8 + +-0.4 + +1.11e-16 + +0.4 + +0.8 + +1.2 + +1.6 + +2 + + + + + + + + + + + + diff --git a/doc/graphs/central_b_spline_second_derivatives.svg b/doc/graphs/central_b_spline_second_derivatives.svg new file mode 100644 index 000000000..7bd209f04 --- /dev/null +++ b/doc/graphs/central_b_spline_second_derivatives.svg @@ -0,0 +1,54 @@ + + + +B-spline Second Derivatives + + + + +-1.615 + +-1.243 + +-0.8699 + +-0.4971 + +-0.1243 + +0.2485 + +0.6213 + +0.9941 + +-1.6 + +-1.2 + +-0.8 + +-0.4 + +1.11e-16 + +0.4 + +0.8 + +1.2 + +1.6 + +2 + + + + + + + + + + diff --git a/doc/sf/cardinal_b_splines.qbk b/doc/sf/cardinal_b_splines.qbk index f4e2f3938..79e61c49a 100644 --- a/doc/sf/cardinal_b_splines.qbk +++ b/doc/sf/cardinal_b_splines.qbk @@ -54,6 +54,8 @@ A basic usage is as follows: [$../graphs/central_b_splines.svg] +[$../graphs/central_b_spline_derivatives.svg] +[$../graphs/central_b_spline_second_derivatives.svg] [h3 Caveats] diff --git a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp index 1347cbc30..684fb76b9 100644 --- a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp @@ -91,7 +91,7 @@ Real b3_spline_double_prime(Real x) if (x < 1) { - return (3*boost::math::constants::half()*x - 2) + x*(3*boost::math::constants::half()); + return 3*x - 2; } if (x < 2) { diff --git a/include/boost/math/special_functions/cardinal_b_spline.hpp b/include/boost/math/special_functions/cardinal_b_spline.hpp index 4e40058ed..5c86d33ff 100644 --- a/include/boost/math/special_functions/cardinal_b_spline.hpp +++ b/include/boost/math/special_functions/cardinal_b_spline.hpp @@ -116,19 +116,26 @@ Real cardinal_b_spline_prime(Real x) } } + if (n==1) + { + if (x==0) + { + return Real(0); + } + if (x==1) + { + return -Real(1)/Real(2); + } + return Real(-1); + } + + Real supp_max = (n+1)/Real(2); if (x >= supp_max) { return Real(0); } - if (n==1) { - if (x==0) { - return Real(0); - } - return Real(-1); - } - // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): std::array v; Real z = x + 1 - supp_max; @@ -154,13 +161,56 @@ Real cardinal_b_spline_prime(Real x) return v[1] - v[0]; } -template -auto forward_cardinal_b_spline(Real x) + +template +Real cardinal_b_spline_double_prime(Real x) { - if constexpr (std::is_integral::value) + static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integer types."); + static_assert(n >= 3); + + if (x < 0) { - return forward_cardinal_b_spline(x); + // All B-splines are even functions, so second derivatives are even: + return cardinal_b_spline_double_prime(-x); } + + + Real supp_max = (n+1)/Real(2); + if (x >= supp_max) + { + return Real(0); + } + + // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): + std::array v; + Real z = x + 1 - supp_max; + for (unsigned i = 0; i < n; ++i) + { + v[i] = detail::B1(z); + z += 1; + } + + Real smx = supp_max - x; + for (unsigned j = 2; j <= n - 2; ++j) + { + Real a = (j + 1 - smx); + Real b = smx; + for(unsigned k = 0; k <= n - j; ++k) + { + v[k] = (a*v[k+1] + b*v[k])/Real(j); + a += 1; + b -= 1; + } + } + + return v[2] - 2*v[1] + v[0]; +} + + +template +Real forward_cardinal_b_spline(Real x) +{ + static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integral types."); return cardinal_b_spline(x - (n+1)/Real(2)); } diff --git a/test/cardinal_b_spline_test.cpp b/test/cardinal_b_spline_test.cpp index 20c5c5dc7..04ad3a0de 100644 --- a/test/cardinal_b_spline_test.cpp +++ b/test/cardinal_b_spline_test.cpp @@ -10,8 +10,8 @@ #include #include #include -#include #include +#include #ifdef BOOST_HAS_FLOAT128 #include using boost::multiprecision::float128; @@ -21,6 +21,8 @@ using std::abs; using boost::math::cardinal_b_spline; using boost::math::cardinal_b_spline_prime; using boost::math::forward_cardinal_b_spline; +using boost::math::cardinal_b_spline_double_prime; + template void test_box() @@ -69,7 +71,15 @@ void test_hat() { std::cerr << " Problem at t = " << t << "\n"; } - if (t < 0) { + if (t == -Real(1)) { + if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0)) { + std::cout << " Problem at t = " << t << "\n"; + } + } + else if (t == Real(1)) { + CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0); + } + else if (t < 0) { CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0); } else if (t == 0) { @@ -105,10 +115,33 @@ void test_quadratic() return (2-t1*t1 -t2*t2)/2; }; + auto b2_prime = [&](Real x)->Real { + Real absx = abs(x); + Real signx = 1; + if (x < 0) { + signx = -1; + } + if (absx >= 3/Real(2)) { + return Real(0); + } + if (absx >= 1/Real(2)) { + return (absx - 3/Real(2))*signx; + } + return -2*absx*signx; + }; + + 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); + expected = b2_prime(t); + + if (!CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<2>(t), 0)) + { + std::cerr << " Problem at t = " << t << "\n"; + } + } } @@ -126,6 +159,18 @@ void test_cubic() expected = Real(0); computed = cardinal_b_spline<3, Real>(2); CHECK_ULP_CLOSE(expected, computed, 0); + + Real h = 1/Real(256); + for (Real t = -4; t <= 4; t += h) { + expected = boost::math::detail::b3_spline_prime(t); + computed = cardinal_b_spline_prime<3>(t); + CHECK_ULP_CLOSE(expected, computed, 0); + expected = boost::math::detail::b3_spline_double_prime(t); + computed = cardinal_b_spline_double_prime<3>(t); + if (!CHECK_ULP_CLOSE(expected, computed, 0)) { + std::cerr << " Problem at t = " << t << "\n"; + } + } } template @@ -149,6 +194,23 @@ void test_quintic() } +template +void test_b_spline_derivatives() +{ + Real h = 1/Real(256); + Real supp = (n+Real(1))/Real(2); + for (Real t = supp - 1; t <= supp+1; t+= h) + { + Real expected = cardinal_b_spline(t+Real(1)/Real(2)) - cardinal_b_spline(t - Real(1)/Real(2)); + Real computed = cardinal_b_spline_prime(t); + CHECK_ULP_CLOSE(expected, computed, 2); + + expected = cardinal_b_spline(t+1) - 2*cardinal_b_spline(t) + cardinal_b_spline(t-1); + computed = cardinal_b_spline_double_prime(t); + CHECK_ULP_CLOSE(expected, computed, 2); + } +} + template void test_partition_of_unity() { @@ -199,6 +261,23 @@ int main() test_partition_of_unity<5, double>(); test_partition_of_unity<6, double>(); + test_b_spline_derivatives<3, double>(); + test_b_spline_derivatives<4, double>(); + test_b_spline_derivatives<5, double>(); + test_b_spline_derivatives<6, double>(); + test_b_spline_derivatives<7, double>(); + test_b_spline_derivatives<8, double>(); + test_b_spline_derivatives<9, double>(); + + test_b_spline_derivatives<3, long double>(); + test_b_spline_derivatives<4, long double>(); + test_b_spline_derivatives<5, long double>(); + test_b_spline_derivatives<6, long double>(); + test_b_spline_derivatives<7, long double>(); + test_b_spline_derivatives<8, long double>(); + test_b_spline_derivatives<9, long double>(); + + #ifdef BOOST_HAS_FLOAT128 test_box(); test_hat(); From edbb25327a2aba95df1d82ba3dbe6f7e2b08426f Mon Sep 17 00:00:00 2001 From: NAThompson Date: Mon, 12 Aug 2019 10:11:58 -0400 Subject: [PATCH 3/4] Cardinal B-spline derivatives: Fix typo in test that made it much less powerful. --- test/cardinal_b_spline_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/cardinal_b_spline_test.cpp b/test/cardinal_b_spline_test.cpp index 04ad3a0de..0ec92cc95 100644 --- a/test/cardinal_b_spline_test.cpp +++ b/test/cardinal_b_spline_test.cpp @@ -199,15 +199,15 @@ void test_b_spline_derivatives() { Real h = 1/Real(256); Real supp = (n+Real(1))/Real(2); - for (Real t = supp - 1; t <= supp+1; t+= h) + for (Real t = -supp - 1; t <= supp+1; t+= h) { Real expected = cardinal_b_spline(t+Real(1)/Real(2)) - cardinal_b_spline(t - Real(1)/Real(2)); Real computed = cardinal_b_spline_prime(t); - CHECK_ULP_CLOSE(expected, computed, 2); + CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits::epsilon()); expected = cardinal_b_spline(t+1) - 2*cardinal_b_spline(t) + cardinal_b_spline(t-1); computed = cardinal_b_spline_double_prime(t); - CHECK_ULP_CLOSE(expected, computed, 2); + CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits::epsilon()); } } From 3291ac891bb99258555369d3cb3e6bbb10628ffd Mon Sep 17 00:00:00 2001 From: NAThompson Date: Tue, 13 Aug 2019 07:54:07 -0400 Subject: [PATCH 4/4] Fix static assert with no message on C++11. --- include/boost/math/special_functions/cardinal_b_spline.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/cardinal_b_spline.hpp b/include/boost/math/special_functions/cardinal_b_spline.hpp index 5c86d33ff..dfb4cf835 100644 --- a/include/boost/math/special_functions/cardinal_b_spline.hpp +++ b/include/boost/math/special_functions/cardinal_b_spline.hpp @@ -166,7 +166,7 @@ template Real cardinal_b_spline_double_prime(Real x) { static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integer types."); - static_assert(n >= 3); + static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required."); if (x < 0) {