Periodic data scheme

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.12.0)
 project(
     BSplineInterpolation
-    VERSION 1.3.0
+    VERSION 1.3.1
     HOMEPAGE_URL "https://github.com/12ff54e/BSplineInterpolation.git"
     LANGUAGES CXX)
 include(CTest)
@@ -42,12 +42,7 @@ install(
     FILES ${PROJECT_BINARY_DIR}/${PROJECT_NAME}${VerPostfix}.cmake
     DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake)
 
-# Set default build type to release.
-if(NOT CMAKE_BUILD_TYPE)
-  set(CMAKE_BUILD_TYPE Release)
-endif()
-
-# Requires c++17 and no compiler extension, but not mandetory.
+# Requires c++11, do not need any compiler extension.
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_EXTENSIONS OFF)
 

README.md

@@ -42,6 +42,18 @@ InterpolationFunction<double, 2> func(
 ```
 Note: this project follows [Semantic Version 2.0.0](https://semver.org/) so the interface will be compatible within one major version.
 
+### Modify Library's Behavior
+
+You can control the bahavior of this library by defining the following macros:
+
+- `INTP_TRACE`, `INTP_DEBUG`: Defines levels of logging, with decreasing verbosity.
+- `INTP_MULTITHREAD`: Multi-threading in solving control points.
+- `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.
+
+By default all of the above macros is not defined.
+
 ### Matlab Interface
 
 Check matlab/Example.m for instructions. Please notice that matlab/bspline.mexw64 is compiled for windows platform, you may need to compile your own MEX file from bspline.cpp (basically a wrapper) on other platforms.
@@ -53,7 +65,7 @@ Check matlab/Example.m for instructions. Please notice that matlab/bspline.mexw6
 
 ## Known Issues
 
-- No
+- Not found yet
 
 ## Future Plan
 

src/include/BSpline.hpp

@@ -6,11 +6,10 @@
 #include <cmath>        // fmod
 #include <functional>   // ref
 #include <iterator>     // distance
-#include <stdexcept>    // range_error
 #include <type_traits>  // is_same, is_arithmatic
 #include <vector>
 
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
 #include <iostream>
 #endif
 
@@ -122,7 +121,7 @@ class BSpline {
                 std::fmod(x - range(dim_ind).first, period) +
                 (x < range(dim_ind).first ? period : knot_type{});
         }
-#ifdef _TRACE
+#ifdef INTP_TRACE
         if ((*iter > x || *(iter + 1) < x) && x >= range(dim_ind).first &&
             x <= range(dim_ind).second) {
             std::cout << "[TRACE] knot hint miss at dim = " << dim_ind
@@ -230,12 +229,11 @@ class BSpline {
               (knot_iter_pairs.second)[-static_cast<int>(order_) - 1])...},
           buf_size_(util::pow(order_ + 1, dim)) {
         for (size_type d = 0; d < dim; ++d) {
-            if (knots_[d].size() - control_points_.dim_size(d) !=
-                (periodicity_[d] ? 2 * order_ + 1 : order_ + 1)) {
-                throw std::range_error(
-                    "Inconsistency between knot number and control point "
-                    "number.");
-            }
+            INTP_ASSERT(knots_[d].size() - control_points_.dim_size(d) ==
+                            (periodicity_[d] ? 2 * order_ + 1 : order_ + 1),
+                        std::string("Inconsistency between knot number and "
+                                    "control point number at dimension ") +
+                            std::to_string(d));
         }
     }
 
@@ -550,9 +548,11 @@ class BSpline {
         return periodicity_[dim_ind];
     }
 
-    inline size_type order() const { return order_; }
+    inline size_type order() const {
+        return order_;
+    }
 
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
     void debug_output() const {
         std::cout << "\n[DEBUG] Control Points (raw data):\n";
 

src/include/BandMatrix.hpp

@@ -48,14 +48,14 @@ class BandMatrix {
      * @return val_type&
      */
     val_type& operator()(size_type i, size_type j) {
-        CUSTOM_ASSERT(j + p_ >= i && i + q_ >= j,
-                      "Given i and j not in main bands.")
+        INTP_ASSERT(j + p_ >= i && i + q_ >= j,
+                    "Given i and j not in main bands.");
         return bands_(j, i + q_ - j);
     }
 
     val_type operator()(size_type i, size_type j) const {
-        CUSTOM_ASSERT(j + p_ >= i && i + q_ >= j,
-                      "Given i and j not in main bands.")
+        INTP_ASSERT(j + p_ >= i && i + q_ >= j,
+                    "Given i and j not in main bands.");
         return bands_(j, i + q_ - j);
     }
 
@@ -120,17 +120,17 @@ class ExtendedBandMatrix : public BandMatrix<T, Alloc> {
     }
 
     val_type& side_bands_val(size_type i, size_type j) {
-        CUSTOM_ASSERT(j >= std::max(n_ - p_, i + q_ + 1) ||
-                          i >= std::max(n_ - q_, j + p_ + 1),
-                      "Given i and j not in side bands.")
+        INTP_ASSERT(j >= std::max(n_ - p_, i + q_ + 1) ||
+                        i >= std::max(n_ - q_, j + p_ + 1),
+                    "Given i and j not in side bands.");
         return j > i + q_ ? right_side_bands_(i, j + p_ - n_)
                           : bottom_side_bands_(j, i + q_ - n_);
     }
 
     val_type side_bands_val(size_type i, size_type j) const {
-        CUSTOM_ASSERT(j >= std::max(n_ - p_, i + q_ + 1) ||
-                          i >= std::max(n_ - q_, j + p_ + 1),
-                      "Given i and j not in side bands.")
+        INTP_ASSERT(j >= std::max(n_ - p_, i + q_ + 1) ||
+                        i >= std::max(n_ - q_, j + p_ + 1),
+                    "Given i and j not in side bands.");
         return j > i + q_ ? right_side_bands_(i, j + p_ - n_)
                           : bottom_side_bands_(j, i + q_ - n_);
     }

src/include/DedicatedThreadPool.hpp

@@ -30,7 +30,7 @@ class DedicatedThreadPool {
      */
     struct shared_working_queue {
        public:
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
         size_t submitted{};
         size_t executed{};
         size_t stolen{};
@@ -40,7 +40,7 @@ class DedicatedThreadPool {
         shared_working_queue() = default;
         void push(task_type&& task) {
             lock_type lk(deq_mutex);
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
             ++submitted;
 #endif
             deq.push_back(std::move(task));
@@ -50,7 +50,7 @@ class DedicatedThreadPool {
             if (deq.empty()) { return false; }
             task = std::move(deq.back());
             deq.pop_back();
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
             ++executed;
 #endif
             return true;
@@ -60,7 +60,7 @@ class DedicatedThreadPool {
             if (deq.empty()) { return false; }
             task = std::move(deq.front());
             deq.pop_front();
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
             ++stolen;
 #endif
             return true;
@@ -95,7 +95,7 @@ class DedicatedThreadPool {
                 threads.emplace_back(&DedicatedThreadPool::thread_loop, this,
                                      i);
             }
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
             std::cout << "[DEBUG] Thread pool initialized with " << t_num
                       << " threads.\n";
 #endif
@@ -114,7 +114,7 @@ class DedicatedThreadPool {
         should_terminate = true;
         cv.notify_all();
 
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
         std::cout << "\n[DEBUG] The thread pool has " << thread_num()
                   << " threads.\n"
                   << "[DEBUG]    Main queue has " << submitted
@@ -142,7 +142,7 @@ class DedicatedThreadPool {
         {
             lock_type lk(main_queue_mutex);
             main_queue.push(std::move(task));
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
             ++submitted;
 #endif
         }
@@ -219,7 +219,7 @@ class DedicatedThreadPool {
         for (size_t i = 0; i < worker_queues.size() - 1; ++i) {
             const auto idx = (thread_idx + i + 1) % worker_queues.size();
             if (worker_queues[idx]->try_steal(task)) {
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
                 ++(worker_queue_ptr->stealing);
 #endif
                 return true;
@@ -254,7 +254,7 @@ class DedicatedThreadPool {
         }
     }
 
-#ifdef _DEBUG
+#ifdef INTP_DEBUG
     size_t submitted{};
 #endif
     bool should_terminate = false;  // Tells threads to stop looking for tasks

src/include/Interpolation.hpp

@@ -4,6 +4,10 @@
 #include <cmath>  // ceil
 #include <initializer_list>
 
+#ifdef INTP_TRACE
+#define INTP_DEBUG
+#endif
+
 #include "InterpolationTemplate.hpp"
 
 namespace intp {
@@ -78,16 +82,21 @@ class InterpolationFunction {  // TODO: Add integration
                         DimArray<typename spline_type::KnotContainer>&,
                         std::pair<T_, T_> x_range) {
         uniform_[dim_ind] = true;
-        const size_type n = mesh_dimension.dim_size(dim_ind);
+        const auto periodic = periodicity(dim_ind);
+        const size_type n = mesh_dimension.dim_size(dim_ind)
+#ifdef INTP_PERIODIC_NO_DUMMY_POINT
+                            + (periodic ? 1 : 0)
+#endif
+            ;
         dx_[dim_ind] =
             (x_range.second - x_range.first) / static_cast<coord_type>(n - 1);
 
         const size_t extra =
-            periodicity(dim_ind) ? 2 * order_ + (1 - order_ % 2) : order_ + 1;
+            periodic ? 2 * order_ + (1 - order_ % 2) : order_ + 1;
 
         std::vector<coord_type> xs(n + extra, x_range.first);
 
-        if (periodicity(dim_ind)) {
+        if (periodic) {
             for (size_type i = 0; i < xs.size(); ++i) {
                 xs[i] = x_range.first +
                         (static_cast<coord_type>(i) -
@@ -120,23 +129,30 @@ class InterpolationFunction {  // TODO: Add integration
         DimArray<typename spline_type::KnotContainer>& input_coords,
         std::pair<T_, T_> x_range) {
         uniform_[dim_ind] = false;
+        const auto periodic = periodicity(dim_ind);
+
         const size_type n{static_cast<size_type>(
             std::distance(x_range.first, x_range.second))};
-        if (n != mesh_dimension.dim_size(dim_ind)) {
-            throw std::range_error(
-                std::string("Inconsistency between knot number and "
-                            "interpolated value number at dimension ") +
-                std::to_string(dim_ind));
-        }
+        INTP_ASSERT(n == mesh_dimension.dim_size(dim_ind)
+#ifdef INTP_PERIODIC_NO_DUMMY_POINT
+                             + (periodic ? 1 : 0)
+#endif
+                        ,
+                    std::string("Inconsistency between knot number and "
+                                "interpolated value number at dimension ") +
+                        std::to_string(dim_ind));
+#ifndef INTP_ENABLE_ASSERTION
+        // suppress unused parameter warning
+        (void)mesh_dimension;
+#endif
         typename spline_type::KnotContainer xs(
-            periodicity(dim_ind) ? n + 2 * order_ + (1 - order_ % 2)
-                                 : n + order_ + 1);
+            periodic ? n + 2 * order_ + (1 - order_ % 2) : n + order_ + 1);
 
         auto& input_coord = input_coords[dim_ind];
         input_coord.reserve(n);
         // The x_range may be given by input iterators, which can not be
         // multi-passed.
-        if (periodicity(dim_ind)) {
+        if (periodic) {
             // In periodic case, the knots are data points, shifted by half of
             // local grid size if spline order is odd.
 
@@ -183,18 +199,15 @@ class InterpolationFunction {  // TODO: Add integration
             // last knot is same as last input coordinate
             input_coord.emplace_back(r_knot);
         }
-#ifdef _DEBUG
         // check whether input coordinates is monotonic
         for (std::size_t i = 0; i < input_coord.size() - 1; ++i) {
-            if (input_coord[i + 1] <= input_coord[i]) {
-                throw std::range_error(
-                    std::string("Given coordinate is not monotonically "
-                                "increasing at dimension ") +
-                    std::to_string(dim_ind));
-            }
+            INTP_ASSERT(input_coord[i + 1] > input_coord[i],
+                        std::string("Given coordinate is not monotonically "
+                                    "increasing at dimension ") +
+                            std::to_string(dim_ind));
         }
-#endif
-#ifdef _TRACE
+
+#ifdef INTP_TRACE
         std::cout << "[TRACE] Nonuniform knots along dimension" << dim_ind
                   << ":\n";
         for (auto& c : xs) { std::cout << "[TRACE] " << c << '\n'; }
@@ -452,7 +465,9 @@ class InterpolationFunction {  // TODO: Add integration
         return spline_.periodicity(dim_ind);
     }
 
-    inline bool uniform(size_type dim_ind) const { return uniform_[dim_ind]; }
+    inline bool uniform(size_type dim_ind) const {
+        return uniform_[dim_ind];
+    }
 
     const std::pair<typename spline_type::knot_type,
                     typename spline_type::knot_type>&
@@ -469,7 +484,9 @@ class InterpolationFunction {  // TODO: Add integration
         return spline_;
     }
 
-    size_type order() const { return order_; }
+    size_type order() const {
+        return order_;
+    }
 };
 
 template <typename T = double, typename U = double>
@@ -485,10 +502,17 @@ class InterpolationFunction1D : public InterpolationFunction<T, size_t{1}, U> {
         : InterpolationFunction1D(
               std::make_pair(typename base::coord_type{},
                              static_cast<typename base::coord_type>(
-                                 (f_range.second - f_range.first) - 1)),
+                                 f_range.second - f_range.first
+#ifdef INTP_PERIODIC_NO_DUMMY_POINT
+                                 - (periodicity ? 0 : 1)
+#else
+                                 - 1
+#endif
+                                     )),
               f_range,
               order_,
-              periodicity) {}
+              periodicity) {
+    }
 
     template <typename C1, typename C2, typename InputIter>
     InterpolationFunction1D(std::pair<C1, C2> x_range,