Weight in cell

.github/workflows/cmake.yml

@@ -33,5 +33,5 @@ jobs:
       working-directory: ${{github.workspace}}/build
       # Execute tests defined by the CMake configuration.  
       # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
-      run: ctest --output-on-failure -R intp-*
+      run: ctest --output-on-failure -L intp-regular
 

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.12.0)
 project(
     BSplineInterpolation
-    VERSION 1.3.1
+    VERSION 1.4.0
     HOMEPAGE_URL "https://github.com/12ff54e/BSplineInterpolation.git"
     LANGUAGES CXX)
 include(CTest)

README.md

@@ -51,6 +51,7 @@ You can control the bahavior of this library by defining the following macros:
 - `INTP_ENABLE_ASSERTION`: Add assertion to ensure data to be interplolated is valid.
 - `INTP_PERIODIC_NO_DUMMY_POINT`: Whether to accept dummy point when interpolating periodic data. If this macro is defined, number of data point is one less than the specified dimension, since the last point is implicitly set as the same of the first one.
 - `INTP_STACK_ALLOCATOR`: Use stack allocator in some small dynamic memory allocation scenario.
+- `INTP_CELL_LAYOUT`: Store weights in cell layout redundantly to speed up evaluation.
 
 By default all of the above macros is not defined.
 

src/include/BSpline.hpp

@@ -9,6 +9,10 @@
 #include <type_traits>  // is_same, is_arithmatic
 #include <vector>
 
+#ifdef INTP_MULTITHREAD
+#include "DedicatedThreadPool.hpp"
+#endif
+
 #ifdef INTP_DEBUG
 #include <iostream>
 #endif
@@ -33,12 +37,21 @@ class BSpline {
     using knot_type = U;
 
     using KnotContainer = std::vector<knot_type>;
+
     using ControlPointContainer =
         Mesh<val_type,
              D,
              util::default_init_allocator<
                  val_type,
                  AlignedAllocator<val_type, Alignment::AVX>>>;
+#ifdef INTP_CELL_LAYOUT
+    using ControlPointCellContainer =
+        Mesh<val_type,
+             D + 1,
+             util::default_init_allocator<
+                 val_type,
+                 AlignedAllocator<val_type, Alignment::AVX>>>;
+#endif
 
     using BaseSpline = std::vector<knot_type>;
     using diff_type = typename KnotContainer::iterator::difference_type;
@@ -156,45 +169,6 @@ class BSpline {
         return {get_knot_iter(indices, coords.first, coords.second)...};
     }
 
-   private:
-    size_type order_;
-
-    DimArray<bool> periodicity_;
-
-    DimArray<KnotContainer> knots_;
-    ControlPointContainer control_points_;
-
-    DimArray<std::pair<knot_type, knot_type>> range_;
-
-    size_type buf_size_;
-
-    // maximum stack buffer size
-    // This buffer is for storing weights when calculating spline derivative
-    // value.
-    constexpr static size_type MAX_BUF_SIZE_ = 1000;
-
-    // auxiliary methods
-
-    /**
-     * @brief Calculate base spline value of each dimension
-     *
-     * @tparam Coords knot_type ...
-     * @tparam indices 0, 1, ...
-     * @param knot_iters an array of knot iters
-     * @param spline_order an array of spline order
-     * @param coords a bunch of coordinates
-     */
-    template <typename... Coords, size_type... indices>
-    inline DimArray<BaseSpline> calc_base_spline_vals(
-        util::index_sequence<indices...>,
-        const DimArray<knot_const_iterator>& knot_iters,
-        const DimArray<size_type>& spline_order,
-        Coords... coords) const {
-        return {base_spline_value(indices, knot_iters[indices], coords,
-                                  spline_order[indices])...};
-    }
-
-   public:
     /**
      * @brief Construct a new BSpline object, with periodicity of each dimension
      * specified.
@@ -217,19 +191,23 @@ class BSpline {
     template <typename... InputIters>
     BSpline(size_type spline_order,
             DimArray<bool> periodicity,
-            ControlPointContainer ctrl_points,
+            ControlPointContainer ctrl_pts,
             std::pair<InputIters, InputIters>... knot_iter_pairs)
         : order_(spline_order),
           periodicity_(periodicity),
           knots_{
               KnotContainer(knot_iter_pairs.first, knot_iter_pairs.second)...},
-          control_points_(std::move(ctrl_points)),
+#ifdef INTP_CELL_LAYOUT
+          control_points_(generate_cell_layout(ctrl_pts)),
+#else
+          control_points_(std::move(ctrl_pts)),
+#endif
           range_{std::make_pair(
               (knot_iter_pairs.first)[order_],
               (knot_iter_pairs.second)[-static_cast<int>(order_) - 1])...},
           buf_size_(util::pow(order_ + 1, dim)) {
         for (size_type d = 0; d < dim; ++d) {
-            INTP_ASSERT(knots_[d].size() - control_points_.dim_size(d) ==
+            INTP_ASSERT(knots_[d].size() - ctrl_pts.dim_size(d) ==
                             (periodicity_[d] ? 2 * order_ + 1 : order_ + 1),
                         std::string("Inconsistency between knot number and "
                                     "control point number at dimension ") +
@@ -258,14 +236,20 @@ class BSpline {
             knots_[dim_ind][knots_[dim_ind].size() - order_ - (2 - order_ % 2)];
     }
 
+#ifdef INTP_CELL_LAYOUT
+    void load_ctrlPts(const ControlPointContainer& control_points) {
+        control_points_ = generate_cell_layout(control_points);
+    }
+#else
     template <typename C>
     typename std::enable_if<
         std::is_same<typename std::remove_reference<C>::type,
                      ControlPointContainer>::value,
         void>::type
-    load_ctrlPts(C&& _control_points) {
-        control_points_ = std::forward<C>(_control_points);
+    load_ctrlPts(C&& control_points) {
+        control_points_ = std::forward<C>(control_points);
     }
+#endif
 
     /**
      * @brief Get spline value at given pairs of coordinate and position hint
@@ -298,12 +282,26 @@ class BSpline {
 
         // combine control points and basic spline values to get spline value
         val_type v{};
+#ifdef INTP_CELL_LAYOUT
+        std::array<size_type, dim + 1> ind_arr{};
+        for (size_type d = 0; d < dim; ++d) {
+            ind_arr[d] = static_cast<size_type>(
+                             distance(knots_begin(d), knot_iters[d])) -
+                         order_;
+        }
+        auto cell_iter = control_points_.begin(dim, ind_arr);
         for (size_type i = 0; i < buf_size_; ++i) {
-            DimArray<size_type> ind_arr;
+            knot_type coef = 1;
             for (size_type d = 0, combined_ind = i; d < dim; ++d) {
-                ind_arr[d] = combined_ind % (order_ + 1);
+                coef *= base_spline_values_1d[d][combined_ind % (order_ + 1)];
                 combined_ind /= (order_ + 1);
             }
+            v += coef * (*cell_iter++);
+        }
+#else
+        MeshDimension<dim> local_mesh_dim(order_ + 1);
+        for (size_type i = 0; i < buf_size_; ++i) {
+            DimArray<size_type> ind_arr = local_mesh_dim.dimwise_indices(i);
 
             knot_type coef = 1;
             for (size_type d = 0; d < dim; ++d) {
@@ -329,6 +327,7 @@ class BSpline {
 
             v += coef * control_points_(ind_arr);
         }
+#endif
 
         return v;
     }
@@ -402,14 +401,39 @@ class BSpline {
 #endif
 
         // get local control points and basic spline values
+
+#ifdef INTP_CELL_LAYOUT
+        std::array<size_type, dim + 1> ind_arr{};
+        for (size_type d = 0; d < dim; ++d) {
+            ind_arr[d] = static_cast<size_type>(
+                             distance(knots_begin(d), knot_iters[d])) -
+                         order_;
+        }
+        auto cell_iter = control_points_.begin(dim, ind_arr);
         for (size_type i = 0; i < buf_size_; ++i) {
             DimArray<size_type> local_ind_arr{};
             for (size_type d = 0, combined_ind = i; d < dim; ++d) {
                 local_ind_arr[d] = combined_ind % (order_ + 1);
                 combined_ind /= (order_ + 1);
             }
 
-            val_type coef = 1;
+            knot_type coef = 1;
+            for (size_type d = 0; d < dim; ++d) {
+                coef *= base_spline_values_1d[d][local_ind_arr[d]];
+            }
+
+            local_spline_val(local_ind_arr) = coef;
+            local_control_points(local_ind_arr) = *cell_iter++;
+        }
+#else
+        for (size_type i = 0; i < buf_size_; ++i) {
+            DimArray<size_type> local_ind_arr{};
+            for (size_type d = 0, combined_ind = i; d < dim; ++d) {
+                local_ind_arr[d] = combined_ind % (order_ + 1);
+                combined_ind /= (order_ + 1);
+            }
+
+            knot_type coef = 1;
             DimArray<size_type> ind_arr{};
             for (size_type d = 0; d < dim; ++d) {
                 coef *= base_spline_values_1d[d][local_ind_arr[d]];
@@ -422,7 +446,8 @@ class BSpline {
                                         knots_begin(d), knot_iters[d])) -
                                         order_);
 
-                // check periodicity, put out-of-right-boundary index to left
+                // check periodicity, put out-of-right-boundary index to
+                // left
                 if (periodicity_[d]) {
                     ind_arr[d] %= control_points_.dim_size(d);
                 }
@@ -431,6 +456,7 @@ class BSpline {
             local_spline_val(local_ind_arr) = coef;
             local_control_points(local_ind_arr) = control_points_(ind_arr);
         }
+#endif
 
         for (size_type d = 0; d < dim; ++d) {
             if (spline_order[d] == order_) { continue; }
@@ -571,6 +597,133 @@ class BSpline {
         std::cout << '\n';
     }
 #endif
+
+   private:
+    size_type order_;
+
+    DimArray<bool> periodicity_;
+
+    DimArray<KnotContainer> knots_;
+#ifdef INTP_CELL_LAYOUT
+    ControlPointCellContainer control_points_;
+#else
+    ControlPointContainer control_points_;
+#endif
+
+    DimArray<std::pair<knot_type, knot_type>> range_;
+
+    size_type buf_size_;
+
+    // maximum stack buffer size
+    // This buffer is for storing weights when calculating spline derivative
+    // value.
+    constexpr static size_type MAX_BUF_SIZE_ = 1000;
+
+    // auxiliary methods
+
+    /**
+     * @brief Calculate base spline value of each dimension
+     *
+     * @tparam Coords knot_type ...
+     * @tparam indices 0, 1, ...
+     * @param knot_iters an array of knot iters
+     * @param spline_order an array of spline order
+     * @param coords a bunch of coordinates
+     */
+    template <typename... Coords, size_type... indices>
+    inline DimArray<BaseSpline> calc_base_spline_vals(
+        util::index_sequence<indices...>,
+        const DimArray<knot_const_iterator>& knot_iters,
+        const DimArray<size_type>& spline_order,
+        Coords... coords) const {
+        return {base_spline_value(indices, knot_iters[indices], coords,
+                                  spline_order[indices])...};
+    }
+
+#ifdef INTP_CELL_LAYOUT
+    ControlPointCellContainer generate_cell_layout(
+        const ControlPointContainer& ctrl_pts) const {
+        MeshDimension<dim + 1> cell_container_dim(
+            util::pow(order_ + 1, dim - 1));
+        for (size_type d = 0; d < dim; ++d) {
+            cell_container_dim.dim_size(d) = ctrl_pts.dim_size(d) -
+                                             (d == dim - 1 ? 0 : order_) +
+                                             (periodicity(d) ? order_ : 0);
+        }
+
+        // Size of the last two dimension of control point cell container. The
+        // (dim-1)th dimention of cell container has the same length (order
+        // points more in periodic case) as the last dimension (which is also
+        // the (dim-1)th dimension) of the origin container, while other
+        // dimensions are #order shorter, except the last dimension with length
+        // (order+1)^(dim-1).
+        const auto line_size = cell_container_dim.dim_size(dim - 1) *
+                               cell_container_dim.dim_size(dim);
+
+        ControlPointCellContainer control_point_cell(cell_container_dim);
+        // size of hyperplane orthogonal to last 2 dimension
+        const size_type hyper_surface_size =
+            control_point_cell.size() / line_size;
+
+        auto fill_cell = [&](size_type begin, size_type end) {
+            // iterate over hyperplane
+            for (size_type h_ind = begin; h_ind < end; ++h_ind) {
+                const auto ind_arr_on_hyper_surface =
+                    control_point_cell.dimension().dimwise_indices(h_ind *
+                                                                   line_size);
+                const auto line_begin =
+                    control_point_cell.begin(dim - 1, ind_arr_on_hyper_surface);
+                const auto line_end =
+                    control_point_cell.end(dim - 1, ind_arr_on_hyper_surface);
+                // iterate along the (dim-1)th dimension
+                for (auto iter = line_begin; iter != line_end; ++iter) {
+                    auto cell_ind_arr = control_point_cell.iter_indices(iter);
+                    MeshDimension<dim - 1> cell_dim(order_ + 1);
+                    // iterate the (dim)th dimension
+                    for (size_type i = 0; i < cell_dim.size(); ++i) {
+                        auto local_shift_ind = cell_dim.dimwise_indices(i);
+                        cell_ind_arr[dim] = i;
+                        auto cp_ind_arr =
+                            typename ControlPointContainer::index_type{};
+                        for (size_type d = 0; d < dim; ++d) {
+                            cp_ind_arr[d] =
+                                (cell_ind_arr[d] +
+                                 (d == dim - 1
+                                      ? 0
+                                      : local_shift_ind[dim - 2 - d])) %
+                                ctrl_pts.dim_size(d);
+                        }
+                        control_point_cell(cell_ind_arr) = ctrl_pts(cp_ind_arr);
+                    }
+                }
+            }
+        };
+#ifdef INTP_MULTITHREAD
+        const size_type block_num = static_cast<size_type>(
+            std::sqrt(static_cast<double>(hyper_surface_size)));
+        const size_type task_per_block = block_num == 0
+                                             ? hyper_surface_size
+                                             : hyper_surface_size / block_num;
+
+        auto& thread_pool = DedicatedThreadPool<void>::get_instance(8);
+        std::vector<std::future<void>> res;
+
+        for (size_type i = 0; i < block_num; ++i) {
+            res.push_back(thread_pool.queue_task([=]() {
+                fill_cell(i * task_per_block, (i + 1) * task_per_block);
+            }));
+        }
+        // main thread deals with the remaining part in case hyper_surface_size
+        // not divisible by thread_num
+        fill_cell(block_num * task_per_block, hyper_surface_size);
+        // wait for all tasks are complete
+        for (auto&& f : res) { f.get(); }
+#else
+        fill_cell(0, hyper_surface_size);
+#endif  // INTP_MULTITHREAD
+        return control_point_cell;
+    }
+#endif  // INTP_CELL_LAYOUT
 };
 
 }  // namespace intp