From d4544ba2184d260d5315d47f668b6734124687ff Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 26 Sep 2023 09:30:27 -0400 Subject: [PATCH 1/9] Compact radial basis function --- ...nterpDetailsCompactRadialBasisFunction.hpp | 94 +++++++++++++++++++ test/CMakeLists.txt | 13 ++- test/tstInterpDetailsCRBF.cpp | 63 +++++++++++++ 3 files changed, 165 insertions(+), 5 deletions(-) create mode 100644 src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp create mode 100644 test/tstInterpDetailsCRBF.cpp diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp new file mode 100644 index 000000000..3049fbcdd --- /dev/null +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -0,0 +1,94 @@ +/**************************************************************************** + * Copyright (c) 2023 by the ArborX authors * + * All rights reserved. * + * * + * This file is part of the ArborX library. ArborX is * + * distributed under a BSD 3-clause license. For the licensing terms see * + * the LICENSE file in the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef ARBORX_INTERP_DETAILS_COMPACT_RADIAL_BASIS_FUNCTION_HPP +#define ARBORX_INTERP_DETAILS_COMPACT_RADIAL_BASIS_FUNCTION_HPP + +#include + +namespace ArborX::Interpolation +{ + +#define CRBF_DECL(NAME) \ + namespace Details \ + { \ + template \ + struct NAME; \ + } \ + \ + namespace CRBF \ + { \ + template \ + static constexpr Details::NAME NAME{}; \ + } + +#define CRBF_DEF(NAME, N, FUNC) \ + namespace Details \ + { \ + template <> \ + struct NAME \ + { \ + template \ + KOKKOS_INLINE_FUNCTION static constexpr T apply(T const y) \ + { \ + T const x = Kokkos::min(Kokkos::abs(y), T(1)); \ + return Kokkos::abs(FUNC); \ + } \ + \ + template \ + KOKKOS_INLINE_FUNCTION constexpr T operator()(T const y) const \ + { \ + return NAME::apply(y); \ + } \ + }; \ + } + +CRBF_DECL(Wendland) +CRBF_DEF(Wendland, 0, (1 - x) * (1 - x)) +CRBF_DEF(Wendland, 2, (1 - x) * (1 - x) * (1 - x) * (1 - x) * (4 * x + 1)) +CRBF_DEF(Wendland, 4, + (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * + (35 * x * x + 18 * x + 3)) +CRBF_DEF(Wendland, 6, + (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * + (1 - x) * (32 * x * x * x + 25 * x * x + 8 * x + 1)) + +CRBF_DECL(Wu) +CRBF_DEF(Wu, 2, + (1 - x) * (1 - x) * (1 - x) * (1 - x) * + (3 * x * x * x + 12 * x + 16 * x + 4)) +CRBF_DEF(Wu, 4, + (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * + (5 * x * x * x * x * x + 30 * x * x * x * x + 72 * x * x * x + + 82 * x * x + 36 * x + 6)) + +CRBF_DECL(Buhmann) +CRBF_DEF(Buhmann, 2, + (x == T(0)) + ? T(1) / 6 + : 2 * x * x * x * x * Kokkos::log(x) - 7 * x * x * x * x / 2 + + 16 * x * x * x / 3 - 2 * x * x + T(1) / 6) +CRBF_DEF(Buhmann, 3, + 1 * x * x * x * x * x * x * x * x - 84 * x * x * x * x * x * x / 5 + + 1024 * x * x * x * x * Kokkos::sqrt(x) / 5 - 378 * x * x * x * x + + 1024 * x * x * x * Kokkos::sqrt(x) / 5 - 84 * x * x / 5 + 1) +CRBF_DEF(Buhmann, 4, + 99 * x * x * x * x * x * x * x * x / 35 - 132 * x * x * x * x * x * x + + 9216 * x * x * x * x * x * Kokkos::sqrt(x) / 35 - + 11264 * x * x * x * x * Kokkos::sqrt(x) / 35 + + 198 * x * x * x * x - 396 * x * x / 35 + 1) + +#undef CRBF_DEF +#undef CRBF_DECL + +} // namespace ArborX::Interpolation + +#endif \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d2a8a3cde..b3fc95687 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -236,11 +236,14 @@ target_link_libraries(ArborX_Test_BoostAdapters.exe PRIVATE ArborX Boost::unit_t target_compile_definitions(ArborX_Test_BoostAdapters.exe PRIVATE BOOST_TEST_DYN_LINK) add_test(NAME ArborX_Test_BoostAdapters COMMAND ArborX_Test_BoostAdapters.exe) -add_executable(ArborX_Test_InterpDetailsSVD.exe tstInterpDetailsSVD.cpp utf_main.cpp) -target_link_libraries(ArborX_Test_InterpDetailsSVD.exe PRIVATE ArborX Boost::unit_test_framework) -target_compile_definitions(ArborX_Test_InterpDetailsSVD.exe PRIVATE BOOST_TEST_DYN_LINK) -target_include_directories(ArborX_Test_InterpDetailsSVD.exe PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) -add_test(NAME ArborX_Test_InterpDetailsSVD COMMAND ArborX_Test_InterpDetailsSVD.exe) +add_executable(ArborX_Test_InterpDetailsMLS.exe + tstInterpDetailsSVD.cpp + tstInterpDetailsCRBF.cpp + utf_main.cpp) +target_link_libraries(ArborX_Test_InterpDetailsMLS.exe PRIVATE ArborX Boost::unit_test_framework) +target_compile_definitions(ArborX_Test_InterpDetailsMLS.exe PRIVATE BOOST_TEST_DYN_LINK) +target_include_directories(ArborX_Test_InterpDetailsMLS.exe PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) +add_test(NAME ArborX_Test_InterpDetailsMLS COMMAND ArborX_Test_InterpDetailsMLS.exe) if(ARBORX_ENABLE_HEADER_SELF_CONTAINMENT_TESTS) add_subdirectory(headers_self_contained) diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCRBF.cpp new file mode 100644 index 000000000..c7f638a34 --- /dev/null +++ b/test/tstInterpDetailsCRBF.cpp @@ -0,0 +1,63 @@ +/**************************************************************************** + * Copyright (c) 2023 by the ArborX authors * + * All rights reserved. * + * * + * This file is part of the ArborX library. ArborX is * + * distributed under a BSD 3-clause license. For the licensing terms see * + * the LICENSE file in the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include + +#include "BoostTest_CUDA_clang_workarounds.hpp" +#include + +#include + +using Functions = std::tuple, + ArborX::Interpolation::Details::Wendland<2>, + ArborX::Interpolation::Details::Wendland<4>, + ArborX::Interpolation::Details::Wendland<6>, + ArborX::Interpolation::Details::Wu<2>, + ArborX::Interpolation::Details::Wu<4>, + ArborX::Interpolation::Details::Buhmann<2>, + ArborX::Interpolation::Details::Buhmann<3>, + ArborX::Interpolation::Details::Buhmann<4>>; + +template +void makeCase(T tol = 1e-5) +{ + static constexpr int range = 15; + + for (int i = 0; i < range; i += 1) + { + T in_unit = T(i) / range; + T out_unit = i + 1; + T zero = 0; + + // ]-inf; -1] + BOOST_TEST(CRBF::apply(-out_unit) >= zero); + BOOST_TEST(CRBF::apply(-out_unit) == zero, + boost::test_tools::tolerance(tol)); + + // ]-1; 1[ + BOOST_TEST(CRBF::apply(-in_unit) >= zero); + BOOST_TEST(CRBF::apply(in_unit) >= zero); + BOOST_TEST(CRBF::apply(in_unit) == CRBF::apply(-in_unit), + boost::test_tools::tolerance(tol)); + + // [1; +inf[ + BOOST_TEST(CRBF::apply(out_unit) >= zero); + BOOST_TEST(CRBF::apply(out_unit) == zero, + boost::test_tools::tolerance(tol)); + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(crbf, Func, Functions) +{ + makeCase(); + makeCase(); + makeCase(); +} \ No newline at end of file From 32396dc40d2b2fd99ba06d5d34be9e5156027aa6 Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 26 Sep 2023 09:44:09 -0400 Subject: [PATCH 2/9] Polynomial basis --- .../ArborX_InterpDetailsPolynomialBasis.hpp | 152 ++++++++++++++++++ test/CMakeLists.txt | 1 + test/tstInterpDetailsPolyBasis.cpp | 72 +++++++++ 3 files changed, 225 insertions(+) create mode 100644 src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp create mode 100644 test/tstInterpDetailsPolyBasis.cpp diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp new file mode 100644 index 000000000..06f1c7256 --- /dev/null +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -0,0 +1,152 @@ +/**************************************************************************** + * Copyright (c) 2023 by the ArborX authors * + * All rights reserved. * + * * + * This file is part of the ArborX library. ArborX is * + * distributed under a BSD 3-clause license. For the licensing terms see * + * the LICENSE file in the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef ARBORX_INTERP_DETAILS_POLYNOMIAL_BASIS_HPP +#define ARBORX_INTERP_DETAILS_POLYNOMIAL_BASIS_HPP + +#include + +#include + +namespace ArborX::Interpolation +{ + +namespace Details +{ + +// Polynomial basis is computed in the same way as the following example +// diagrams. Each cell is the product of the variable times everything on the +// left cell and above. +// +// For example, the "degree 3 / variable y" has size 3 because it multiplies +// every element of "degree 2 / variable y" and "degree 2 / variable x". +// +// Quadratic | 3D || Cubic | 2D +// ------------+----- || ------------------+----- +// degree | var || degree | var +// ---+---+----+ || ---+---+----+-----+ +// 0 | 1 | 2 | || 0 | 1 | 2 | 3 | +// ---+---+----+ || ---+---+----+-----+ +// ---+---+----+ || ---+---+----+-----+ +// 1 | | | || 1 | | | | +// ---+---+----+----- || ---+---+----+-----+----- +// | x | | x || | x | | | x +// | | xx | || | | xx | | +// +---+----+----- || | | | xxx | +// | y | | y || +---+----+-----+----- +// | | xy | || | y | | | y +// | | yy | || | | xy | | +// +---+----+----- || | | yy | | +// | z | | z || | | | xxy | +// | | xz | || | | | xyy | +// | | yz | || | | | yyy | +// | | zz | || + +// The size of each cell is 1 if it is at degree 1 or at variable x +// And the sum of the cells' size to the left and above. (1 alone is not +// counted). This essentially creates Pascal's triangle where line n corresponds +// to the "dim + deg = n"-th diagonal +template +KOKKOS_FUNCTION constexpr auto polynomialBasisCellSizes() +{ + static_assert(Deg != 0 && Dim != 0, + "Unable to compute cell sizes for a constant polynomial basis"); + + struct + { + std::size_t arr[Deg][Dim]{}; + } result; + auto &arr = result.arr; + + for (std::size_t dim = 0; dim < Dim; dim++) + arr[0][dim] = 1; + + for (std::size_t deg = 0; deg < Deg; deg++) + arr[deg][0] = 1; + + for (std::size_t deg = 1; deg < Deg; deg++) + for (std::size_t dim = 1; dim < Dim; dim++) + arr[deg][dim] = arr[deg - 1][dim] + arr[deg][dim - 1]; + + return result; +} + +// This returns the size of the polynomial basis, which is the sum of all the +// cells' sizes and 1. +template +KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() +{ + if constexpr (Deg != 0 && Dim != 0) + { + auto [arr] = polynomialBasisCellSizes(); + std::size_t size = 1; + + for (std::size_t deg = 0; deg < Deg; deg++) + for (std::size_t dim = 0; dim < Dim; dim++) + size += arr[deg][dim]; + + return size; + } + else + { + return 1; + } +} + +// This builds the array as described above +template +KOKKOS_FUNCTION auto polynomialBasis(Point const &p) +{ + using value_t = std::decay_t; + Kokkos::Array()> arr{}; + arr[0] = value_t(1); + + if constexpr (Deg != 0 && Dim != 0) + { + // Cannot use struct binding with constexpr + static constexpr auto cell_sizes_struct = + polynomialBasisCellSizes(); + static constexpr auto &cell_sizes = cell_sizes_struct.arr; + + std::size_t prev_col = 0; + std::size_t curr_col = 1; + + for (std::size_t deg = 0; deg < Deg; deg++) + { + std::size_t loc_offset = curr_col; + for (std::size_t dim = 0; dim < Dim; dim++) + { + // copy the previous column and multply by p[dim] + for (std::size_t i = 0; i < cell_sizes[deg][dim]; i++) + arr[loc_offset + i] = arr[prev_col + i] * p[dim]; + + loc_offset += cell_sizes[deg][dim]; + } + + prev_col = curr_col; + curr_col = loc_offset; + } + } + + return arr; +} + +template +struct PolynomialDegree : std::integral_constant +{}; + +} // namespace Details +template +static constexpr Details::PolynomialDegree PolynomialDegree{}; + +} // namespace ArborX::Interpolation + +#endif \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b3fc95687..3ff9ebd9c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -239,6 +239,7 @@ add_test(NAME ArborX_Test_BoostAdapters COMMAND ArborX_Test_BoostAdapters.exe) add_executable(ArborX_Test_InterpDetailsMLS.exe tstInterpDetailsSVD.cpp tstInterpDetailsCRBF.cpp + tstInterpDetailsPolyBasis.cpp utf_main.cpp) target_link_libraries(ArborX_Test_InterpDetailsMLS.exe PRIVATE ArborX Boost::unit_test_framework) target_compile_definitions(ArborX_Test_InterpDetailsMLS.exe PRIVATE BOOST_TEST_DYN_LINK) diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp new file mode 100644 index 000000000..3ef418593 --- /dev/null +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -0,0 +1,72 @@ +/**************************************************************************** + * Copyright (c) 2023 by the ArborX authors * + * All rights reserved. * + * * + * This file is part of the ArborX library. ArborX is * + * distributed under a BSD 3-clause license. For the licensing terms see * + * the LICENSE file in the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include "ArborX_EnableViewComparison.hpp" +#include + +#include + +BOOST_AUTO_TEST_CASE(polynomial_basis_column_sizes) +{ + using view = Kokkos::View; + + auto [arr0] = + ArborX::Interpolation::Details::polynomialBasisCellSizes<5, 3>(); + std::size_t ref0[3][5] = { + {1, 1, 1, 1, 1}, {1, 2, 3, 4, 5}, {1, 3, 6, 10, 15}}; + view arr0_view(&arr0[0][0], 3, 5); + view ref0_view(&ref0[0][0], 3, 5); + ARBORX_MDVIEW_TEST(arr0_view, ref0_view); + + auto [arr1] = + ArborX::Interpolation::Details::polynomialBasisCellSizes<2, 3>(); + std::size_t ref1[3][2] = {{1, 1}, {1, 2}, {1, 3}}; + view arr1_view(&arr1[0][0], 3, 2); + view ref1_view(&ref1[0][0], 3, 2); + ARBORX_MDVIEW_TEST(arr1_view, ref1_view); +} + +BOOST_AUTO_TEST_CASE(polynomial_basis_size) +{ + auto sz0 = ArborX::Interpolation::Details::polynomialBasisSize<5, 3>(); + BOOST_TEST(sz0 == 56); + + auto sz1 = ArborX::Interpolation::Details::polynomialBasisSize<2, 3>(); + BOOST_TEST(sz1 == 10); + + auto sz2 = ArborX::Interpolation::Details::polynomialBasisSize<3, 0>(); + BOOST_TEST(sz2 == 1); + + auto sz3 = ArborX::Interpolation::Details::polynomialBasisSize<0, 3>(); + BOOST_TEST(sz3 == 1); +} + +BOOST_AUTO_TEST_CASE(polynomial_basis) +{ + using view = Kokkos::View; + + double point0[5] = {1, 2, 3, 4, 5}; + auto arr0 = ArborX::Interpolation::Details::polynomialBasis<5, 3>(point0); + double ref0[56] = {1, 1, 2, 3, 4, 5, 1, 2, 4, 3, 6, 9, 4, 8, + 12, 16, 5, 10, 15, 20, 25, 1, 2, 4, 8, 3, 6, 12, + 9, 18, 27, 4, 8, 16, 12, 24, 36, 16, 32, 48, 64, 5, + 10, 20, 15, 30, 45, 20, 40, 60, 80, 25, 50, 75, 100, 125}; + view arr0_view(arr0.data(), 56); + view ref0_view(&ref0[0], 56); + ARBORX_MDVIEW_TEST(arr0_view, ref0_view); + + double point1[2] = {-2, 0}; + auto arr1 = ArborX::Interpolation::Details::polynomialBasis<2, 3>(point1); + double ref1[10] = {1, -2, 0, 4, 0, 0, -8, 0, 0, 0}; + view arr1_view(arr1.data(), 10); + view ref1_view(&ref1[0], 10); + ARBORX_MDVIEW_TEST(arr1_view, ref1_view); +} \ No newline at end of file From 6ffe87cc6488db6d47789a23000a95c342a7219f Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Thu, 28 Sep 2023 10:30:06 -0400 Subject: [PATCH 3/9] Fixes and extra comments --- ...nterpDetailsCompactRadialBasisFunction.hpp | 28 +++------------- .../ArborX_InterpDetailsPolynomialBasis.hpp | 26 +++++++++++---- test/tstInterpDetailsCRBF.cpp | 32 +++++++++---------- test/tstInterpDetailsPolyBasis.cpp | 6 ++-- 4 files changed, 45 insertions(+), 47 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index 3049fbcdd..622998cf0 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -14,42 +14,24 @@ #include -namespace ArborX::Interpolation +namespace ArborX::Interpolation::CRBF { #define CRBF_DECL(NAME) \ - namespace Details \ - { \ template \ - struct NAME; \ - } \ - \ - namespace CRBF \ - { \ - template \ - static constexpr Details::NAME NAME{}; \ - } + struct NAME; #define CRBF_DEF(NAME, N, FUNC) \ - namespace Details \ - { \ template <> \ struct NAME \ { \ template \ - KOKKOS_INLINE_FUNCTION static constexpr T apply(T const y) \ + KOKKOS_INLINE_FUNCTION static constexpr T evaluate(T const y) \ { \ T const x = Kokkos::min(Kokkos::abs(y), T(1)); \ return Kokkos::abs(FUNC); \ } \ - \ - template \ - KOKKOS_INLINE_FUNCTION constexpr T operator()(T const y) const \ - { \ - return NAME::apply(y); \ - } \ - }; \ - } + }; CRBF_DECL(Wendland) CRBF_DEF(Wendland, 0, (1 - x) * (1 - x)) @@ -89,6 +71,6 @@ CRBF_DEF(Buhmann, 4, #undef CRBF_DEF #undef CRBF_DECL -} // namespace ArborX::Interpolation +} // namespace ArborX::Interpolation::CRBF #endif \ No newline at end of file diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index 06f1c7256..fa3dcd80a 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -49,11 +49,24 @@ namespace Details // | | xz | || | | | xyy | // | | yz | || | | | yyy | // | | zz | || +// +// This generates: || This generates: +// [1, x, y, z, xx, || [1, x, y, xx, xy, yy, +// xy, yy, xz, yz, || xxx, xxy, xyy, yyy] +// zz] || // The size of each cell is 1 if it is at degree 1 or at variable x // And the sum of the cells' size to the left and above. (1 alone is not // counted). This essentially creates Pascal's triangle where line n corresponds // to the "dim + deg = n"-th diagonal +// +// Given the previous two diagrams, the two 2D arrays would be: +// +// Deg \ Dim | x | y | z || Deg \ Dim | x | y +// ----------+---+---+--- || ----------+---+--- +// 1 | 1 | 1 | 1 || 1 | 1 | 1 +// 2 | 1 | 2 | 3 || 2 | 1 | 2 +// || 3 | 1 | 3 template KOKKOS_FUNCTION constexpr auto polynomialBasisCellSizes() { @@ -81,6 +94,9 @@ KOKKOS_FUNCTION constexpr auto polynomialBasisCellSizes() // This returns the size of the polynomial basis, which is the sum of all the // cells' sizes and 1. +// +// Given the previous two diagrams and 2D arrays, both the Quadratic/3D and +// Cubic/2D would have 10 elements in total. template KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() { @@ -103,7 +119,7 @@ KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() // This builds the array as described above template -KOKKOS_FUNCTION auto polynomialBasis(Point const &p) +KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) { using value_t = std::decay_t; Kokkos::Array()> arr{}; @@ -111,7 +127,7 @@ KOKKOS_FUNCTION auto polynomialBasis(Point const &p) if constexpr (Deg != 0 && Dim != 0) { - // Cannot use struct binding with constexpr + // Cannot use structured binding with constexpr static constexpr auto cell_sizes_struct = polynomialBasisCellSizes(); static constexpr auto &cell_sizes = cell_sizes_struct.arr; @@ -139,14 +155,12 @@ KOKKOS_FUNCTION auto polynomialBasis(Point const &p) return arr; } +} // namespace Details + template struct PolynomialDegree : std::integral_constant {}; -} // namespace Details -template -static constexpr Details::PolynomialDegree PolynomialDegree{}; - } // namespace ArborX::Interpolation #endif \ No newline at end of file diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCRBF.cpp index c7f638a34..af1b2bd24 100644 --- a/test/tstInterpDetailsCRBF.cpp +++ b/test/tstInterpDetailsCRBF.cpp @@ -16,15 +16,15 @@ #include -using Functions = std::tuple, - ArborX::Interpolation::Details::Wendland<2>, - ArborX::Interpolation::Details::Wendland<4>, - ArborX::Interpolation::Details::Wendland<6>, - ArborX::Interpolation::Details::Wu<2>, - ArborX::Interpolation::Details::Wu<4>, - ArborX::Interpolation::Details::Buhmann<2>, - ArborX::Interpolation::Details::Buhmann<3>, - ArborX::Interpolation::Details::Buhmann<4>>; +using Functions = std::tuple, + ArborX::Interpolation::CRBF::Wendland<2>, + ArborX::Interpolation::CRBF::Wendland<4>, + ArborX::Interpolation::CRBF::Wendland<6>, + ArborX::Interpolation::CRBF::Wu<2>, + ArborX::Interpolation::CRBF::Wu<4>, + ArborX::Interpolation::CRBF::Buhmann<2>, + ArborX::Interpolation::CRBF::Buhmann<3>, + ArborX::Interpolation::CRBF::Buhmann<4>>; template void makeCase(T tol = 1e-5) @@ -38,19 +38,19 @@ void makeCase(T tol = 1e-5) T zero = 0; // ]-inf; -1] - BOOST_TEST(CRBF::apply(-out_unit) >= zero); - BOOST_TEST(CRBF::apply(-out_unit) == zero, + BOOST_TEST(CRBF::evaluate(-out_unit) >= zero); + BOOST_TEST(CRBF::evaluate(-out_unit) == zero, boost::test_tools::tolerance(tol)); // ]-1; 1[ - BOOST_TEST(CRBF::apply(-in_unit) >= zero); - BOOST_TEST(CRBF::apply(in_unit) >= zero); - BOOST_TEST(CRBF::apply(in_unit) == CRBF::apply(-in_unit), + BOOST_TEST(CRBF::evaluate(-in_unit) >= zero); + BOOST_TEST(CRBF::evaluate(in_unit) >= zero); + BOOST_TEST(CRBF::evaluate(in_unit) == CRBF::evaluate(-in_unit), boost::test_tools::tolerance(tol)); // [1; +inf[ - BOOST_TEST(CRBF::apply(out_unit) >= zero); - BOOST_TEST(CRBF::apply(out_unit) == zero, + BOOST_TEST(CRBF::evaluate(out_unit) >= zero); + BOOST_TEST(CRBF::evaluate(out_unit) == zero, boost::test_tools::tolerance(tol)); } } diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp index 3ef418593..fd4ad7add 100644 --- a/test/tstInterpDetailsPolyBasis.cpp +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -54,7 +54,8 @@ BOOST_AUTO_TEST_CASE(polynomial_basis) using view = Kokkos::View; double point0[5] = {1, 2, 3, 4, 5}; - auto arr0 = ArborX::Interpolation::Details::polynomialBasis<5, 3>(point0); + auto arr0 = + ArborX::Interpolation::Details::evaluatePolynomialBasis<5, 3>(point0); double ref0[56] = {1, 1, 2, 3, 4, 5, 1, 2, 4, 3, 6, 9, 4, 8, 12, 16, 5, 10, 15, 20, 25, 1, 2, 4, 8, 3, 6, 12, 9, 18, 27, 4, 8, 16, 12, 24, 36, 16, 32, 48, 64, 5, @@ -64,7 +65,8 @@ BOOST_AUTO_TEST_CASE(polynomial_basis) ARBORX_MDVIEW_TEST(arr0_view, ref0_view); double point1[2] = {-2, 0}; - auto arr1 = ArborX::Interpolation::Details::polynomialBasis<2, 3>(point1); + auto arr1 = + ArborX::Interpolation::Details::evaluatePolynomialBasis<2, 3>(point1); double ref1[10] = {1, -2, 0, 4, 0, 0, -8, 0, 0, 0}; view arr1_view(arr1.data(), 10); view ref1_view(&ref1[0], 10); From 96af593e50b78a613e1a88472381e0b301c92701 Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Mon, 2 Oct 2023 13:21:51 -0400 Subject: [PATCH 4/9] Better explanations and crbf test --- .../ArborX_InterpDetailsPolynomialBasis.hpp | 96 +++++++++---------- test/tstInterpDetailsCRBF.cpp | 4 +- test/tstInterpDetailsPolyBasis.cpp | 6 +- 3 files changed, 50 insertions(+), 56 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index fa3dcd80a..5f25cb6fb 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -22,56 +22,52 @@ namespace ArborX::Interpolation namespace Details { -// Polynomial basis is computed in the same way as the following example -// diagrams. Each cell is the product of the variable times everything on the -// left cell and above. +// The goal of these functions is to evaluate the polynomial basis at any degree +// and dimension. For example, the polynomial basis of degree 2 evaluated at a +// point {x, y, z} would be [1, x, y, z, xx, xy, yy, xz, yz, zz]. // -// For example, the "degree 3 / variable y" has size 3 because it multiplies -// every element of "degree 2 / variable y" and "degree 2 / variable x". +// To compute it, the list of values is sliced against its degree and highest +// dimension. From the previous list, we would have the following sublists: +// - [1] at degree 0 +// - [x], [y] and [z] at degree 1 +// - [xx], [xy, yy] and [xz, yz, zz] at degree 2 +// One can then infer a recursive pattern that will be used to build the list: +// - [1] at degree 0 +// - x*[1], y*[1] and z*[1] at degree 1 +// - x*[x], y*[x, y] and z*[x, y, z] at degree 2 +// So, given a slice at degree n and dimension u, its values would be the +// product of all the slices of degree n-1 and of dimension u or less. // -// Quadratic | 3D || Cubic | 2D -// ------------+----- || ------------------+----- -// degree | var || degree | var -// ---+---+----+ || ---+---+----+-----+ -// 0 | 1 | 2 | || 0 | 1 | 2 | 3 | -// ---+---+----+ || ---+---+----+-----+ -// ---+---+----+ || ---+---+----+-----+ -// 1 | | | || 1 | | | | -// ---+---+----+----- || ---+---+----+-----+----- -// | x | | x || | x | | | x -// | | xx | || | | xx | | -// +---+----+----- || | | | xxx | -// | y | | y || +---+----+-----+----- -// | | xy | || | y | | | y -// | | yy | || | | xy | | -// +---+----+----- || | | yy | | -// | z | | z || | | | xxy | -// | | xz | || | | | xyy | -// | | yz | || | | | yyy | -// | | zz | || -// -// This generates: || This generates: -// [1, x, y, z, xx, || [1, x, y, xx, xy, yy, -// xy, yy, xz, yz, || xxx, xxy, xyy, yyy] -// zz] || - -// The size of each cell is 1 if it is at degree 1 or at variable x -// And the sum of the cells' size to the left and above. (1 alone is not -// counted). This essentially creates Pascal's triangle where line n corresponds -// to the "dim + deg = n"-th diagonal -// -// Given the previous two diagrams, the two 2D arrays would be: +// As another example, if we can take the polynomial basis of degree 3 evaluated +// at a point {x, y} would be [1, x, y, xx, xy, yy, xxx, xxy, xyy, yyy]. Its +// slices are: +// - [1] at degree 0 +// - [x] and [y] at degree 1 +// - [xx] and [xy, yy] at degree 2 +// - [xxx] and [xxy, xyy, yyy] at degree 3 +// And its recursive pattern: +// - [1] at degree 0 +// - x*[1] and y*[1] at degree 1 +// - x*[x] and y*[x, y] at degree 2 +// - x*[xx] and y*[xx, xy, yy] at degree 3 + +// This function returns the lengths of the slices for a given maximum degree +// and dimension. +// Given a slice at degree n and dimension u, its length would be the sum of all +// the slices' lengths of degree n-1 and of same or smaller dimension. // +// Given the previous two examples, the two arrays would be: // Deg \ Dim | x | y | z || Deg \ Dim | x | y // ----------+---+---+--- || ----------+---+--- // 1 | 1 | 1 | 1 || 1 | 1 | 1 // 2 | 1 | 2 | 3 || 2 | 1 | 2 // || 3 | 1 | 3 template -KOKKOS_FUNCTION constexpr auto polynomialBasisCellSizes() +KOKKOS_FUNCTION constexpr auto polynomialBasisSliceLengths() { - static_assert(Deg != 0 && Dim != 0, - "Unable to compute cell sizes for a constant polynomial basis"); + static_assert( + Deg != 0 && Dim != 0, + "Unable to compute slice lengths for a constant polynomial basis"); struct { @@ -93,16 +89,14 @@ KOKKOS_FUNCTION constexpr auto polynomialBasisCellSizes() } // This returns the size of the polynomial basis, which is the sum of all the -// cells' sizes and 1. -// -// Given the previous two diagrams and 2D arrays, both the Quadratic/3D and -// Cubic/2D would have 10 elements in total. +// slices lengths and 1. +// Both of the previous examples will result in a polynomial basis of size 10 template KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() { if constexpr (Deg != 0 && Dim != 0) { - auto [arr] = polynomialBasisCellSizes(); + auto [arr] = polynomialBasisSliceLengths(); std::size_t size = 1; for (std::size_t deg = 0; deg < Deg; deg++) @@ -117,7 +111,7 @@ KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() } } -// This builds the array as described above +// This creates the list by building each slices in-place template KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) { @@ -128,9 +122,9 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) if constexpr (Deg != 0 && Dim != 0) { // Cannot use structured binding with constexpr - static constexpr auto cell_sizes_struct = - polynomialBasisCellSizes(); - static constexpr auto &cell_sizes = cell_sizes_struct.arr; + static constexpr auto slice_lengths_struct = + polynomialBasisSliceLengths(); + static constexpr auto &slice_lengths = slice_lengths_struct.arr; std::size_t prev_col = 0; std::size_t curr_col = 1; @@ -141,10 +135,10 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) for (std::size_t dim = 0; dim < Dim; dim++) { // copy the previous column and multply by p[dim] - for (std::size_t i = 0; i < cell_sizes[deg][dim]; i++) + for (std::size_t i = 0; i < slice_lengths[deg][dim]; i++) arr[loc_offset + i] = arr[prev_col + i] * p[dim]; - loc_offset += cell_sizes[deg][dim]; + loc_offset += slice_lengths[deg][dim]; } prev_col = curr_col; diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCRBF.cpp index af1b2bd24..5f1a6ab7c 100644 --- a/test/tstInterpDetailsCRBF.cpp +++ b/test/tstInterpDetailsCRBF.cpp @@ -43,8 +43,8 @@ void makeCase(T tol = 1e-5) boost::test_tools::tolerance(tol)); // ]-1; 1[ - BOOST_TEST(CRBF::evaluate(-in_unit) >= zero); - BOOST_TEST(CRBF::evaluate(in_unit) >= zero); + BOOST_TEST(CRBF::evaluate(-in_unit) > zero); + BOOST_TEST(CRBF::evaluate(in_unit) > zero); BOOST_TEST(CRBF::evaluate(in_unit) == CRBF::evaluate(-in_unit), boost::test_tools::tolerance(tol)); diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp index fd4ad7add..878ebb18f 100644 --- a/test/tstInterpDetailsPolyBasis.cpp +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -14,12 +14,12 @@ #include -BOOST_AUTO_TEST_CASE(polynomial_basis_column_sizes) +BOOST_AUTO_TEST_CASE(polynomial_basis_slice_lengths) { using view = Kokkos::View; auto [arr0] = - ArborX::Interpolation::Details::polynomialBasisCellSizes<5, 3>(); + ArborX::Interpolation::Details::polynomialBasisSliceLengths<5, 3>(); std::size_t ref0[3][5] = { {1, 1, 1, 1, 1}, {1, 2, 3, 4, 5}, {1, 3, 6, 10, 15}}; view arr0_view(&arr0[0][0], 3, 5); @@ -27,7 +27,7 @@ BOOST_AUTO_TEST_CASE(polynomial_basis_column_sizes) ARBORX_MDVIEW_TEST(arr0_view, ref0_view); auto [arr1] = - ArborX::Interpolation::Details::polynomialBasisCellSizes<2, 3>(); + ArborX::Interpolation::Details::polynomialBasisSliceLengths<2, 3>(); std::size_t ref1[3][2] = {{1, 1}, {1, 2}, {1, 3}}; view arr1_view(&arr1[0][0], 3, 2); view ref1_view(&ref1[0][0], 3, 2); From 6b1c23492c2f03f9f92319289608afcfd7e330e5 Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Mon, 2 Oct 2023 15:33:31 -0400 Subject: [PATCH 5/9] CRBF readability and poly basis comments --- ...nterpDetailsCompactRadialBasisFunction.hpp | 68 +++++++++++-------- .../ArborX_InterpDetailsPolynomialBasis.hpp | 37 ++++++---- 2 files changed, 65 insertions(+), 40 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index 622998cf0..7d9889679 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -14,7 +14,24 @@ #include -namespace ArborX::Interpolation::CRBF +namespace ArborX::Interpolation +{ + +namespace Details +{ + +template +KOKKOS_INLINE_FUNCTION T evaluatePolynomial(T const x, T const (&coeffs)[N]) +{ + T eval = 0; + for (std::size_t i = 0; i < N; i++) + eval = x * eval + coeffs[i]; + return eval; +} + +} // namespace Details + +namespace CRBF { #define CRBF_DECL(NAME) \ @@ -33,44 +50,39 @@ namespace ArborX::Interpolation::CRBF } \ }; +#define CRBF_POLY(...) Details::evaluatePolynomial(x, {__VA_ARGS__}) +#define CRBF_POW(X, N) Kokkos::pow(X, N) + CRBF_DECL(Wendland) -CRBF_DEF(Wendland, 0, (1 - x) * (1 - x)) -CRBF_DEF(Wendland, 2, (1 - x) * (1 - x) * (1 - x) * (1 - x) * (4 * x + 1)) -CRBF_DEF(Wendland, 4, - (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * - (35 * x * x + 18 * x + 3)) -CRBF_DEF(Wendland, 6, - (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * - (1 - x) * (32 * x * x * x + 25 * x * x + 8 * x + 1)) +CRBF_DEF(Wendland, 0, CRBF_POW(1 - x, 2)) +CRBF_DEF(Wendland, 2, CRBF_POW(1 - x, 4) * CRBF_POLY(4, 1)) +CRBF_DEF(Wendland, 4, CRBF_POW(1 - x, 6) * CRBF_POLY(35, 18, 3)) +CRBF_DEF(Wendland, 6, CRBF_POW(1 - x, 6) * CRBF_POLY(32, 25, 8, 1)) CRBF_DECL(Wu) -CRBF_DEF(Wu, 2, - (1 - x) * (1 - x) * (1 - x) * (1 - x) * - (3 * x * x * x + 12 * x + 16 * x + 4)) -CRBF_DEF(Wu, 4, - (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * (1 - x) * - (5 * x * x * x * x * x + 30 * x * x * x * x + 72 * x * x * x + - 82 * x * x + 36 * x + 6)) +CRBF_DEF(Wu, 2, CRBF_POW(1 - x, 4) * CRBF_POLY(3, 12, 16, 4)) +CRBF_DEF(Wu, 4, CRBF_POW(1 - x, 6) * CRBF_POLY(5, 30, 72, 82, 36, 6)) CRBF_DECL(Buhmann) CRBF_DEF(Buhmann, 2, - (x == T(0)) - ? T(1) / 6 - : 2 * x * x * x * x * Kokkos::log(x) - 7 * x * x * x * x / 2 + - 16 * x * x * x / 3 - 2 * x * x + T(1) / 6) + (x == T(0)) ? T(1) / 6 + : CRBF_POLY(12 * Kokkos::log(x) - 21, 32, -12, 0, 1) / 6) CRBF_DEF(Buhmann, 3, - 1 * x * x * x * x * x * x * x * x - 84 * x * x * x * x * x * x / 5 + - 1024 * x * x * x * x * Kokkos::sqrt(x) / 5 - 378 * x * x * x * x + - 1024 * x * x * x * Kokkos::sqrt(x) / 5 - 84 * x * x / 5 + 1) + CRBF_POLY(5, 0, -84, 0, 1024 * Kokkos::sqrt(x) - 1890, + 1024 * Kokkos::sqrt(x), -84, 0, 5) / + 5) CRBF_DEF(Buhmann, 4, - 99 * x * x * x * x * x * x * x * x / 35 - 132 * x * x * x * x * x * x + - 9216 * x * x * x * x * x * Kokkos::sqrt(x) / 35 - - 11264 * x * x * x * x * Kokkos::sqrt(x) / 35 + - 198 * x * x * x * x - 396 * x * x / 35 + 1) + CRBF_POLY(99, 0, -4620, 9216 * Kokkos::sqrt(x), + -11264 * Kokkos::sqrt(x) + 6930, 0, -396, 0, 35) / + 35) +#undef CRBF_POW +#undef CRBF_POLY #undef CRBF_DEF #undef CRBF_DECL -} // namespace ArborX::Interpolation::CRBF +} // namespace CRBF + +} // namespace ArborX::Interpolation #endif \ No newline at end of file diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index 5f25cb6fb..ce789ee9f 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -16,10 +16,7 @@ #include -namespace ArborX::Interpolation -{ - -namespace Details +namespace ArborX::Interpolation::Details { // The goal of these functions is to evaluate the polynomial basis at any degree @@ -36,7 +33,8 @@ namespace Details // - x*[1], y*[1] and z*[1] at degree 1 // - x*[x], y*[x, y] and z*[x, y, z] at degree 2 // So, given a slice at degree n and dimension u, its values would be the -// product of all the slices of degree n-1 and of dimension u or less. +// product of all the slices of degree n-1 and of dimension u or less with the +// coordinate of dimension u. // // As another example, if we can take the polynomial basis of degree 3 evaluated // at a point {x, y} would be [1, x, y, xx, xy, yy, xxx, xxy, xyy, yyy]. Its @@ -50,6 +48,27 @@ namespace Details // - x*[1] and y*[1] at degree 1 // - x*[x] and y*[x, y] at degree 2 // - x*[xx] and y*[xx, xy, yy] at degree 3 +// +// These examples can be represented in 2D as the following tables: +// Quadratic | 3D || Cubic | 2D +// ------------+----- || ------------------+----- +// degree | dim || degree | dim +// ---+---+----+ || ---+---+----+-----+ +// 0 | 1 | 2 | || 0 | 1 | 2 | 3 | +// ===o===o====o || ===o===o====o=====o +// 1 | | || 1 | | +// ---+---+----+----- || ---+---+----+-----+----- +// | x | | x || | x | | | x +// | | xx | || | | xx | | +// +---+----+----- || | | | xxx | +// | y | | y || +---+----+-----+----- +// | | xy | || | y | | | y +// | | yy | || | | xy | | +// +---+----+----- || | | yy | | +// | z | | z || | | | xxy | +// | | xz | || | | | xyy | +// | | yz | || | | | yyy | +// | | zz | || // This function returns the lengths of the slices for a given maximum degree // and dimension. @@ -149,12 +168,6 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) return arr; } -} // namespace Details - -template -struct PolynomialDegree : std::integral_constant -{}; - -} // namespace ArborX::Interpolation +} // namespace ArborX::Interpolation::Details #endif \ No newline at end of file From 4a4b56ac112977d54efecfbd104605e977576bcb Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 3 Oct 2023 14:42:09 -0400 Subject: [PATCH 6/9] Testing CRBF against true function (if available) --- ...nterpDetailsCompactRadialBasisFunction.hpp | 23 ++-- test/ArborX_EnableViewComparison.hpp | 2 +- test/tstInterpDetailsCRBF.cpp | 107 +++++++++++------- 3 files changed, 81 insertions(+), 51 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index 7d9889679..1ad4acc79 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -14,18 +14,21 @@ #include +#include + namespace ArborX::Interpolation { namespace Details { -template -KOKKOS_INLINE_FUNCTION T evaluatePolynomial(T const x, T const (&coeffs)[N]) +template +KOKKOS_INLINE_FUNCTION T +evaluatePolynomial(T const x, std::initializer_list const coeffs) { T eval = 0; - for (std::size_t i = 0; i < N; i++) - eval = x * eval + coeffs[i]; + for (auto const coeff : coeffs) + eval = x * eval + coeff; return eval; } @@ -54,14 +57,14 @@ namespace CRBF #define CRBF_POW(X, N) Kokkos::pow(X, N) CRBF_DECL(Wendland) -CRBF_DEF(Wendland, 0, CRBF_POW(1 - x, 2)) -CRBF_DEF(Wendland, 2, CRBF_POW(1 - x, 4) * CRBF_POLY(4, 1)) -CRBF_DEF(Wendland, 4, CRBF_POW(1 - x, 6) * CRBF_POLY(35, 18, 3)) -CRBF_DEF(Wendland, 6, CRBF_POW(1 - x, 6) * CRBF_POLY(32, 25, 8, 1)) +CRBF_DEF(Wendland, 0, CRBF_POW(CRBF_POLY(-1, 1), 2)) +CRBF_DEF(Wendland, 2, CRBF_POW(CRBF_POLY(-1, 1), 4) * CRBF_POLY(4, 1)) +CRBF_DEF(Wendland, 4, CRBF_POW(CRBF_POLY(-1, 1), 6) * CRBF_POLY(35, 18, 3)) +CRBF_DEF(Wendland, 6, CRBF_POW(CRBF_POLY(-1, 1), 8) * CRBF_POLY(32, 25, 8, 1)) CRBF_DECL(Wu) -CRBF_DEF(Wu, 2, CRBF_POW(1 - x, 4) * CRBF_POLY(3, 12, 16, 4)) -CRBF_DEF(Wu, 4, CRBF_POW(1 - x, 6) * CRBF_POLY(5, 30, 72, 82, 36, 6)) +CRBF_DEF(Wu, 2, CRBF_POW(CRBF_POLY(-1, 1), 4) * CRBF_POLY(3, 12, 16, 4)) +CRBF_DEF(Wu, 4, CRBF_POW(CRBF_POLY(-1, 1), 6) * CRBF_POLY(5, 30, 72, 82, 36, 6)) CRBF_DECL(Buhmann) CRBF_DEF(Buhmann, 2, diff --git a/test/ArborX_EnableViewComparison.hpp b/test/ArborX_EnableViewComparison.hpp index 26b741298..720b440a8 100644 --- a/test/ArborX_EnableViewComparison.hpp +++ b/test/ArborX_EnableViewComparison.hpp @@ -84,7 +84,7 @@ void arborxViewCheck(U const &u, V const &v, std::string const &u_name, } #define ARBORX_MDVIEW_TEST_TOL(VIEWA, VIEWB, TOL) \ - [](decltype(VIEWA) const &u, decltype(VIEWB) const &v) { \ + [&](decltype(VIEWA) const &u, decltype(VIEWB) const &v) { \ auto view_a = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, u); \ auto view_b = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, v); \ \ diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCRBF.cpp index 5f1a6ab7c..7d7f29a1c 100644 --- a/test/tstInterpDetailsCRBF.cpp +++ b/test/tstInterpDetailsCRBF.cpp @@ -9,55 +9,82 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ +#include "ArborX_EnableDeviceTypes.hpp" +#include "ArborX_EnableViewComparison.hpp" #include #include "BoostTest_CUDA_clang_workarounds.hpp" +#include #include -#include - -using Functions = std::tuple, - ArborX::Interpolation::CRBF::Wendland<2>, - ArborX::Interpolation::CRBF::Wendland<4>, - ArborX::Interpolation::CRBF::Wendland<6>, - ArborX::Interpolation::CRBF::Wu<2>, - ArborX::Interpolation::CRBF::Wu<4>, - ArborX::Interpolation::CRBF::Buhmann<2>, - ArborX::Interpolation::CRBF::Buhmann<3>, - ArborX::Interpolation::CRBF::Buhmann<4>>; - -template -void makeCase(T tol = 1e-5) +template +void makeCase(ES const &es, CRBF, TrueFunc const &tf, T tol = 1e-5) { + using View = Kokkos::View; + using HostView = typename View::HostMirror; static constexpr int range = 15; - for (int i = 0; i < range; i += 1) + HostView input("Testing::input", 4 * range); + for (int i = 0; i < range; i++) { - T in_unit = T(i) / range; - T out_unit = i + 1; - T zero = 0; - - // ]-inf; -1] - BOOST_TEST(CRBF::evaluate(-out_unit) >= zero); - BOOST_TEST(CRBF::evaluate(-out_unit) == zero, - boost::test_tools::tolerance(tol)); - - // ]-1; 1[ - BOOST_TEST(CRBF::evaluate(-in_unit) > zero); - BOOST_TEST(CRBF::evaluate(in_unit) > zero); - BOOST_TEST(CRBF::evaluate(in_unit) == CRBF::evaluate(-in_unit), - boost::test_tools::tolerance(tol)); - - // [1; +inf[ - BOOST_TEST(CRBF::evaluate(out_unit) >= zero); - BOOST_TEST(CRBF::evaluate(out_unit) == zero, - boost::test_tools::tolerance(tol)); + input(4 * i + 0) = -i - 1; + input(4 * i + 1) = -T(i) / range; + input(4 * i + 2) = T(i) / range; + input(4 * i + 3) = i + 1; + } + + View eval("Testing::eval", 4 * range); + Kokkos::deep_copy(es, eval, input); + Kokkos::parallel_for( + "Testing::eval_crbf", Kokkos::RangePolicy(es, 0, 4 * range), + KOKKOS_LAMBDA(int const i) { eval(i) = CRBF::evaluate(eval(i)); }); + + if constexpr (!std::is_same_v) + { + HostView reference("Testing::reference", 4 * range); + for (int i = 0; i < range; i++) + { + reference(4 * i + 0) = 0; + reference(4 * i + 1) = tf(T(i) / range); + reference(4 * i + 2) = tf(T(i) / range); + reference(4 * i + 3) = 0; + } + + ARBORX_MDVIEW_TEST_TOL(eval, reference, tol); } + + auto heval = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, eval); + for (int i = 0; i < 4 * range; i++) + BOOST_TEST(heval(i) >= T(0)); } -BOOST_AUTO_TEST_CASE_TEMPLATE(crbf, Func, Functions) -{ - makeCase(); - makeCase(); - makeCase(); -} \ No newline at end of file +#define MAKE_TEST(F, I, TF) \ + BOOST_AUTO_TEST_CASE_TEMPLATE(crbf_##F##_##I, DeviceType, \ + ARBORX_DEVICE_TYPES) \ + { \ + makeCase(typename DeviceType::execution_space{}, \ + ArborX::Interpolation::CRBF::F{}, TF); \ + } + +#define MAKE_TEST_POLY(F, I, POLY) \ + MAKE_TEST(F, I, [](auto x) -> decltype(x) { \ + using boost::math::tools::pow; \ + using poly = boost::math::tools::polynomial; \ + return (POLY)(x); \ + }) + +#define MAKE_TEST_NONE(F, I) MAKE_TEST(F, I, nullptr) + +MAKE_TEST_POLY(Wendland, 0, (pow(poly{1, -1}, 2))) +MAKE_TEST_POLY(Wendland, 2, (pow(poly{1, -1}, 4) * poly{1, 4})) +MAKE_TEST_POLY(Wendland, 4, (pow(poly{1, -1}, 6) * poly{3, 18, 35})) +MAKE_TEST_POLY(Wendland, 6, (pow(poly{1, -1}, 8) * poly{1, 8, 25, 32})) +MAKE_TEST_POLY(Wu, 2, (pow(poly{1, -1}, 4) * poly{4, 16, 12, 3})) +MAKE_TEST_POLY(Wu, 4, (pow(poly{1, -1}, 6) * poly{6, 36, 82, 72, 30, 5})) +MAKE_TEST_NONE(Buhmann, 2) +MAKE_TEST_NONE(Buhmann, 3) +MAKE_TEST_NONE(Buhmann, 4) + +#undef MAKE_TEST_NONE +#undef MAKE_TEST_POLY +#undef MAKE_TEST \ No newline at end of file From e896eefa77d005f8dbd3b9500ddd4d2e5ec04ed8 Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 3 Oct 2023 15:44:23 -0400 Subject: [PATCH 7/9] Comments and renaming --- ...nterpDetailsCompactRadialBasisFunction.hpp | 6 + .../ArborX_InterpDetailsPolynomialBasis.hpp | 109 +++++++----------- test/tstInterpDetailsPolyBasis.cpp | 17 +-- 3 files changed, 58 insertions(+), 74 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index 1ad4acc79..e124b0df7 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -22,6 +22,8 @@ namespace ArborX::Interpolation namespace Details { +// Polynomials are represented with the highest coefficient first. For example, +// x^2 - 2x + 1 would be {1, -2, 1} template KOKKOS_INLINE_FUNCTION T evaluatePolynomial(T const x, std::initializer_list const coeffs) @@ -48,6 +50,10 @@ namespace CRBF template \ KOKKOS_INLINE_FUNCTION static constexpr T evaluate(T const y) \ { \ + /* We force the input to be between 0 and 1. \ + Because CRBF(-a) = CRBF(a) = CRBF(|a|), we take the absolute value \ + and clamp the range to [0, 1] before entering in the definition of \ + the CRBF. */ \ T const x = Kokkos::min(Kokkos::abs(y), T(1)); \ return Kokkos::abs(FUNC); \ } \ diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index ce789ee9f..5d9f00341 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -12,9 +12,9 @@ #ifndef ARBORX_INTERP_DETAILS_POLYNOMIAL_BASIS_HPP #define ARBORX_INTERP_DETAILS_POLYNOMIAL_BASIS_HPP -#include +#include -#include +#include namespace ArborX::Interpolation::Details { @@ -49,109 +49,86 @@ namespace ArborX::Interpolation::Details // - x*[x] and y*[x, y] at degree 2 // - x*[xx] and y*[xx, xy, yy] at degree 3 // -// These examples can be represented in 2D as the following tables: -// Quadratic | 3D || Cubic | 2D -// ------------+----- || ------------------+----- -// degree | dim || degree | dim -// ---+---+----+ || ---+---+----+-----+ -// 0 | 1 | 2 | || 0 | 1 | 2 | 3 | -// ===o===o====o || ===o===o====o=====o -// 1 | | || 1 | | -// ---+---+----+----- || ---+---+----+-----+----- -// | x | | x || | x | | | x -// | | xx | || | | xx | | -// +---+----+----- || | | | xxx | -// | y | | y || +---+----+-----+----- -// | | xy | || | y | | | y -// | | yy | || | | xy | | -// +---+----+----- || | | yy | | -// | z | | z || | | | xxy | -// | | xz | || | | | xyy | -// | | yz | || | | | yyy | -// | | zz | || - -// This function returns the lengths of the slices for a given maximum degree -// and dimension. -// Given a slice at degree n and dimension u, its length would be the sum of all -// the slices' lengths of degree n-1 and of same or smaller dimension. -// -// Given the previous two examples, the two arrays would be: -// Deg \ Dim | x | y | z || Deg \ Dim | x | y -// ----------+---+---+--- || ----------+---+--- -// 1 | 1 | 1 | 1 || 1 | 1 | 1 -// 2 | 1 | 2 | 3 || 2 | 1 | 2 -// || 3 | 1 | 3 -template +// The lengths for the slices would be +// Deg \ Dim | x | y | z +// ----------+---+---+--- +// 1 | 1 | 1 | 1 +// 2 | 1 | 2 | 3 +// 3 | 1 | 3 | 6 +template KOKKOS_FUNCTION constexpr auto polynomialBasisSliceLengths() { static_assert( - Deg != 0 && Dim != 0, + Degree > 0 && DIM > 0, "Unable to compute slice lengths for a constant polynomial basis"); struct { - std::size_t arr[Deg][Dim]{}; + std::size_t arr[Degree][DIM]{}; } result; auto &arr = result.arr; - for (std::size_t dim = 0; dim < Dim; dim++) + for (std::size_t dim = 0; dim < DIM; dim++) arr[0][dim] = 1; - for (std::size_t deg = 0; deg < Deg; deg++) + for (std::size_t deg = 0; deg < Degree; deg++) arr[deg][0] = 1; - for (std::size_t deg = 1; deg < Deg; deg++) - for (std::size_t dim = 1; dim < Dim; dim++) + for (std::size_t deg = 1; deg < Degree; deg++) + for (std::size_t dim = 1; dim < DIM; dim++) arr[deg][dim] = arr[deg - 1][dim] + arr[deg][dim - 1]; return result; } -// This returns the size of the polynomial basis, which is the sum of all the -// slices lengths and 1. -// Both of the previous examples will result in a polynomial basis of size 10 -template -KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() +// This returns the size of the polynomial basis. Counting the constant 1, the +// size would be "Deg + Dim choose Dim" or "Deg + Dim choose Deg". +KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize(std::size_t dim, + std::size_t deg) { - if constexpr (Deg != 0 && Dim != 0) - { - auto [arr] = polynomialBasisSliceLengths(); - std::size_t size = 1; + std::size_t dim_fact = 1; + std::size_t deg_fact = 1; + std::size_t dim_deg_fact = 1; - for (std::size_t deg = 0; deg < Deg; deg++) - for (std::size_t dim = 0; dim < Dim; dim++) - size += arr[deg][dim]; + for (std::size_t i = 2; i <= dim; i++) + dim_fact *= i; - return size; - } - else - { - return 1; - } + for (std::size_t i = 2; i <= deg; i++) + deg_fact *= i; + + for (std::size_t i = 2; i <= dim + deg; i++) + dim_deg_fact *= i; + + return dim_deg_fact / (dim_fact * deg_fact); } // This creates the list by building each slices in-place -template +template KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) { - using value_t = std::decay_t; - Kokkos::Array()> arr{}; + GeometryTraits::check_valid_geometry_traits(Point{}); + static_assert(GeometryTraits::is_point::value, + "point must be a point"); + static constexpr std::size_t DIM = GeometryTraits::dimension_v; + using value_t = typename GeometryTraits::coordinate_type::type; + + Kokkos::Array arr{}; arr[0] = value_t(1); - if constexpr (Deg != 0 && Dim != 0) + if constexpr (Degree > 0 && DIM > 0) { // Cannot use structured binding with constexpr static constexpr auto slice_lengths_struct = - polynomialBasisSliceLengths(); + polynomialBasisSliceLengths(); static constexpr auto &slice_lengths = slice_lengths_struct.arr; std::size_t prev_col = 0; std::size_t curr_col = 1; - for (std::size_t deg = 0; deg < Deg; deg++) + for (std::size_t deg = 0; deg < Degree; deg++) { std::size_t loc_offset = curr_col; - for (std::size_t dim = 0; dim < Dim; dim++) + for (std::size_t dim = 0; dim < DIM; dim++) { // copy the previous column and multply by p[dim] for (std::size_t i = 0; i < slice_lengths[deg][dim]; i++) diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp index 878ebb18f..1b3409c78 100644 --- a/test/tstInterpDetailsPolyBasis.cpp +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -10,6 +10,7 @@ ****************************************************************************/ #include "ArborX_EnableViewComparison.hpp" +#include #include #include @@ -36,16 +37,16 @@ BOOST_AUTO_TEST_CASE(polynomial_basis_slice_lengths) BOOST_AUTO_TEST_CASE(polynomial_basis_size) { - auto sz0 = ArborX::Interpolation::Details::polynomialBasisSize<5, 3>(); + auto sz0 = ArborX::Interpolation::Details::polynomialBasisSize(5, 3); BOOST_TEST(sz0 == 56); - auto sz1 = ArborX::Interpolation::Details::polynomialBasisSize<2, 3>(); + auto sz1 = ArborX::Interpolation::Details::polynomialBasisSize(2, 3); BOOST_TEST(sz1 == 10); - auto sz2 = ArborX::Interpolation::Details::polynomialBasisSize<3, 0>(); + auto sz2 = ArborX::Interpolation::Details::polynomialBasisSize(3, 0); BOOST_TEST(sz2 == 1); - auto sz3 = ArborX::Interpolation::Details::polynomialBasisSize<0, 3>(); + auto sz3 = ArborX::Interpolation::Details::polynomialBasisSize(0, 3); BOOST_TEST(sz3 == 1); } @@ -53,9 +54,9 @@ BOOST_AUTO_TEST_CASE(polynomial_basis) { using view = Kokkos::View; - double point0[5] = {1, 2, 3, 4, 5}; + ArborX::ExperimentalHyperGeometry::Point<5, double> point0 = {1, 2, 3, 4, 5}; auto arr0 = - ArborX::Interpolation::Details::evaluatePolynomialBasis<5, 3>(point0); + ArborX::Interpolation::Details::evaluatePolynomialBasis<3>(point0); double ref0[56] = {1, 1, 2, 3, 4, 5, 1, 2, 4, 3, 6, 9, 4, 8, 12, 16, 5, 10, 15, 20, 25, 1, 2, 4, 8, 3, 6, 12, 9, 18, 27, 4, 8, 16, 12, 24, 36, 16, 32, 48, 64, 5, @@ -64,9 +65,9 @@ BOOST_AUTO_TEST_CASE(polynomial_basis) view ref0_view(&ref0[0], 56); ARBORX_MDVIEW_TEST(arr0_view, ref0_view); - double point1[2] = {-2, 0}; + ArborX::ExperimentalHyperGeometry::Point<2, double> point1 = {-2, 0}; auto arr1 = - ArborX::Interpolation::Details::evaluatePolynomialBasis<2, 3>(point1); + ArborX::Interpolation::Details::evaluatePolynomialBasis<3>(point1); double ref1[10] = {1, -2, 0, 4, 0, 0, -8, 0, 0, 0}; view arr1_view(arr1.data(), 10); view ref1_view(&ref1[0], 10); From 3d0ab88cdabaf2c449e6f787fd6787905b8db052 Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 3 Oct 2023 16:50:16 -0400 Subject: [PATCH 8/9] nvcc fix and dimension assertion --- ...nterpDetailsCompactRadialBasisFunction.hpp | 2 +- .../ArborX_InterpDetailsPolynomialBasis.hpp | 20 +++++++++------ test/tstInterpDetailsCRBF.cpp | 25 ++++++++++++------- test/tstInterpDetailsPolyBasis.cpp | 17 +++++-------- 4 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index e124b0df7..23d12bd14 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -23,7 +23,7 @@ namespace Details { // Polynomials are represented with the highest coefficient first. For example, -// x^2 - 2x + 1 would be {1, -2, 1} +// 3x^2 - 2x + 1 would be {3, -2, 1} template KOKKOS_INLINE_FUNCTION T evaluatePolynomial(T const x, std::initializer_list const coeffs) diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index 5d9f00341..a5cdc8a7f 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -58,8 +58,9 @@ namespace ArborX::Interpolation::Details template KOKKOS_FUNCTION constexpr auto polynomialBasisSliceLengths() { + static_assert(DIM > 0, "Polynomial basis with no dimension is invalid"); static_assert( - Degree > 0 && DIM > 0, + Degree > 0, "Unable to compute slice lengths for a constant polynomial basis"); struct @@ -83,20 +84,22 @@ KOKKOS_FUNCTION constexpr auto polynomialBasisSliceLengths() // This returns the size of the polynomial basis. Counting the constant 1, the // size would be "Deg + Dim choose Dim" or "Deg + Dim choose Deg". -KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize(std::size_t dim, - std::size_t deg) +template +KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() { + static_assert(DIM > 0, "Polynomial basis with no dimension is invalid"); + std::size_t dim_fact = 1; std::size_t deg_fact = 1; std::size_t dim_deg_fact = 1; - for (std::size_t i = 2; i <= dim; i++) + for (std::size_t i = 2; i <= DIM; i++) dim_fact *= i; - for (std::size_t i = 2; i <= deg; i++) + for (std::size_t i = 2; i <= Degree; i++) deg_fact *= i; - for (std::size_t i = 2; i <= dim + deg; i++) + for (std::size_t i = 2; i <= DIM + Degree; i++) dim_deg_fact *= i; return dim_deg_fact / (dim_fact * deg_fact); @@ -111,11 +114,12 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) "point must be a point"); static constexpr std::size_t DIM = GeometryTraits::dimension_v; using value_t = typename GeometryTraits::coordinate_type::type; + static_assert(DIM > 0, "Polynomial basis with no dimension is invalid"); - Kokkos::Array arr{}; + Kokkos::Array()> arr{}; arr[0] = value_t(1); - if constexpr (Degree > 0 && DIM > 0) + if constexpr (Degree > 0) { // Cannot use structured binding with constexpr static constexpr auto slice_lengths_struct = diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCRBF.cpp index 7d7f29a1c..808b98d9c 100644 --- a/test/tstInterpDetailsCRBF.cpp +++ b/test/tstInterpDetailsCRBF.cpp @@ -17,8 +17,8 @@ #include #include -template -void makeCase(ES const &es, CRBF, TrueFunc const &tf, T tol = 1e-5) +template +void makeCase(ES const &es, std::function const &tf, T tol = 1e-5) { using View = Kokkos::View; using HostView = typename View::HostMirror; @@ -39,7 +39,7 @@ void makeCase(ES const &es, CRBF, TrueFunc const &tf, T tol = 1e-5) "Testing::eval_crbf", Kokkos::RangePolicy(es, 0, 4 * range), KOKKOS_LAMBDA(int const i) { eval(i) = CRBF::evaluate(eval(i)); }); - if constexpr (!std::is_same_v) + if (bool(tf)) { HostView reference("Testing::reference", 4 * range); for (int i = 0; i < range; i++) @@ -62,18 +62,25 @@ void makeCase(ES const &es, CRBF, TrueFunc const &tf, T tol = 1e-5) BOOST_AUTO_TEST_CASE_TEMPLATE(crbf_##F##_##I, DeviceType, \ ARBORX_DEVICE_TYPES) \ { \ - makeCase(typename DeviceType::execution_space{}, \ - ArborX::Interpolation::CRBF::F{}, TF); \ + makeCase>( \ + typename DeviceType::execution_space{}, TF); \ } #define MAKE_TEST_POLY(F, I, POLY) \ - MAKE_TEST(F, I, [](auto x) -> decltype(x) { \ + template \ + T func##F##I(T const x) \ + { \ using boost::math::tools::pow; \ - using poly = boost::math::tools::polynomial; \ + using poly = boost::math::tools::polynomial; \ return (POLY)(x); \ - }) + } \ + \ + MAKE_TEST(F, I, func##F##I) + +template +static std::function emptyFunc = {}; -#define MAKE_TEST_NONE(F, I) MAKE_TEST(F, I, nullptr) +#define MAKE_TEST_NONE(F, I) MAKE_TEST(F, I, emptyFunc) MAKE_TEST_POLY(Wendland, 0, (pow(poly{1, -1}, 2))) MAKE_TEST_POLY(Wendland, 2, (pow(poly{1, -1}, 4) * poly{1, 4})) diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp index 1b3409c78..6ae2aaba0 100644 --- a/test/tstInterpDetailsPolyBasis.cpp +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -37,17 +37,12 @@ BOOST_AUTO_TEST_CASE(polynomial_basis_slice_lengths) BOOST_AUTO_TEST_CASE(polynomial_basis_size) { - auto sz0 = ArborX::Interpolation::Details::polynomialBasisSize(5, 3); - BOOST_TEST(sz0 == 56); - - auto sz1 = ArborX::Interpolation::Details::polynomialBasisSize(2, 3); - BOOST_TEST(sz1 == 10); - - auto sz2 = ArborX::Interpolation::Details::polynomialBasisSize(3, 0); - BOOST_TEST(sz2 == 1); - - auto sz3 = ArborX::Interpolation::Details::polynomialBasisSize(0, 3); - BOOST_TEST(sz3 == 1); + BOOST_TEST( + (ArborX::Interpolation::Details::polynomialBasisSize<5, 3>() == 56)); + BOOST_TEST( + (ArborX::Interpolation::Details::polynomialBasisSize<2, 3>() == 10)); + BOOST_TEST( + (ArborX::Interpolation::Details::polynomialBasisSize<3, 0>() == 1)); } BOOST_AUTO_TEST_CASE(polynomial_basis) From 79947bdffb1f5b185ded5bdf3eda3f0e5bec342f Mon Sep 17 00:00:00 2001 From: Yohann Bosqued Date: Tue, 10 Oct 2023 09:44:19 -0400 Subject: [PATCH 9/9] Better basis size computation --- .../ArborX_InterpDetailsPolynomialBasis.hpp | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index a5cdc8a7f..e39a9bc28 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -89,20 +89,10 @@ KOKKOS_FUNCTION constexpr std::size_t polynomialBasisSize() { static_assert(DIM > 0, "Polynomial basis with no dimension is invalid"); - std::size_t dim_fact = 1; - std::size_t deg_fact = 1; - std::size_t dim_deg_fact = 1; - - for (std::size_t i = 2; i <= DIM; i++) - dim_fact *= i; - - for (std::size_t i = 2; i <= Degree; i++) - deg_fact *= i; - - for (std::size_t i = 2; i <= DIM + Degree; i++) - dim_deg_fact *= i; - - return dim_deg_fact / (dim_fact * deg_fact); + std::size_t result = 1; + for (std::size_t k = 0; k < Kokkos::min(DIM, Degree); ++k) + result = result * (DIM + Degree - k) / (k + 1); + return result; } // This creates the list by building each slices in-place