Compact radial basis functions and generic polynomial basis / Utility for moving least squares

src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp

@@ -0,0 +1,97 @@
+/****************************************************************************
+ * 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 <Kokkos_Core.hpp>
+
+#include <initializer_list>
+
+namespace ArborX::Interpolation
+{
+
+namespace Details
+{
+
+// Polynomials are represented with the highest coefficient first. For example,
+// 3x^2 - 2x + 1 would be {3, -2, 1}
+template <typename T>
+KOKKOS_INLINE_FUNCTION T
+evaluatePolynomial(T const x, std::initializer_list<T> const coeffs)
+{
+  T eval = 0;
+  for (auto const coeff : coeffs)
+    eval = x * eval + coeff;
+  return eval;
+}
+
+} // namespace Details
+
+namespace CRBF
+{
+
+#define CRBF_DECL(NAME)                                                        \
+  template <std::size_t>                                                       \
+  struct NAME;
+
+#define CRBF_DEF(NAME, N, FUNC)                                                \
+  template <>                                                                  \
+  struct NAME<N>                                                               \
+  {                                                                            \
+    template <typename T>                                                      \
+    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);                                                \
+    }                                                                          \
+  };
+
+#define CRBF_POLY(...) Details::evaluatePolynomial<T>(x, {__VA_ARGS__})
+#define CRBF_POW(X, N) Kokkos::pow(X, N)
+
+CRBF_DECL(Wendland)
+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(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,
+         (x == T(0)) ? T(1) / 6
+                     : CRBF_POLY(12 * Kokkos::log(x) - 21, 32, -12, 0, 1) / 6)
+CRBF_DEF(Buhmann, 3,
+         CRBF_POLY(5, 0, -84, 0, 1024 * Kokkos::sqrt(x) - 1890,
+                   1024 * Kokkos::sqrt(x), -84, 0, 5) /
+             5)
+CRBF_DEF(Buhmann, 4,
+         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 CRBF
+
+} // namespace ArborX::Interpolation
+
+#endif

src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp

@@ -0,0 +1,154 @@
+/****************************************************************************
+ * 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 <ArborX_GeometryTraits.hpp>
+
+#include <Kokkos_Core.hpp>
+
+namespace ArborX::Interpolation::Details
+{
+
+// 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].
+//
+// 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 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
+// 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
+//
+// 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 <std::size_t DIM, std::size_t Degree>
+KOKKOS_FUNCTION constexpr auto polynomialBasisSliceLengths()
+{
+  static_assert(DIM > 0, "Polynomial basis with no dimension is invalid");
+  static_assert(
+      Degree > 0,
+      "Unable to compute slice lengths for a constant polynomial basis");
+
+  struct
+  {
+    std::size_t arr[Degree][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 < Degree; deg++)
+    arr[deg][0] = 1;
+
+  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. Counting the constant 1, the
+// size would be "Deg + Dim choose Dim" or "Deg + Dim choose Deg".
+template <std::size_t DIM, std::size_t Degree>
+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);
+}
+
+// This creates the list by building each slices in-place
+template <std::size_t Degree, typename Point>
+KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p)
+{
+  GeometryTraits::check_valid_geometry_traits(Point{});
+  static_assert(GeometryTraits::is_point<Point>::value,
+                "point must be a point");
+  static constexpr std::size_t DIM = GeometryTraits::dimension_v<Point>;
+  using value_t = typename GeometryTraits::coordinate_type<Point>::type;
+  static_assert(DIM > 0, "Polynomial basis with no dimension is invalid");
+
+  Kokkos::Array<value_t, polynomialBasisSize<DIM, Degree>()> arr{};
+  arr[0] = value_t(1);
+
+  if constexpr (Degree > 0)
+  {
+    // Cannot use structured binding with constexpr
+    static constexpr auto slice_lengths_struct =
+        polynomialBasisSliceLengths<DIM, Degree>();
+    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 < Degree; 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 < slice_lengths[deg][dim]; i++)
+          arr[loc_offset + i] = arr[prev_col + i] * p[dim];
+
+        loc_offset += slice_lengths[deg][dim];
+      }
+
+      prev_col = curr_col;
+      curr_col = loc_offset;
+    }
+  }
+
+  return arr;
+}
+
+} // namespace ArborX::Interpolation::Details
+
+#endif

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); \
                                                                                \

test/CMakeLists.txt

@@ -236,11 +236,15 @@ 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
+  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)
+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)

test/tstInterpDetailsCRBF.cpp

@@ -0,0 +1,97 @@
+/****************************************************************************
+ * 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_EnableDeviceTypes.hpp"
+#include "ArborX_EnableViewComparison.hpp"
+#include <interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp>
+
+#include "BoostTest_CUDA_clang_workarounds.hpp"
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/test/unit_test.hpp>
+
+template <typename T, typename CRBF, typename ES>
+void makeCase(ES const &es, std::function<T(T const)> const &tf, T tol = 1e-5)
+{
+  using View = Kokkos::View<T *, typename ES::memory_space>;
+  using HostView = typename View::HostMirror;
+  static constexpr int range = 15;
+
+  HostView input("Testing::input", 4 * range);
+  for (int i = 0; i < range; i++)
+  {
+    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>(es, 0, 4 * range),
+      KOKKOS_LAMBDA(int const i) { eval(i) = CRBF::evaluate(eval(i)); });
+
+  if (bool(tf))
+  {
+    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));
+}
+
+#define MAKE_TEST(F, I, TF)                                                    \
+  BOOST_AUTO_TEST_CASE_TEMPLATE(crbf_##F##_##I, DeviceType,                    \
+                                ARBORX_DEVICE_TYPES)                           \
+  {                                                                            \
+    makeCase<double, ArborX::Interpolation::CRBF::F<I>>(                       \
+        typename DeviceType::execution_space{}, TF<double>);                   \
+  }
+
+#define MAKE_TEST_POLY(F, I, POLY)                                             \
+  template <typename T>                                                        \
+  T func##F##I(T const x)                                                      \
+  {                                                                            \
+    using boost::math::tools::pow;                                             \
+    using poly = boost::math::tools::polynomial<T>;                            \
+    return (POLY)(x);                                                          \
+  }                                                                            \
+                                                                               \
+  MAKE_TEST(F, I, func##F##I)
+
+template <typename T>
+static std::function<T(T const)> emptyFunc = {};
+
+#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}))
+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