From ffe82952cc6841f780c288919aa55c9cce7b5a45 Mon Sep 17 00:00:00 2001 From: qzhong <12ff54e@gmail.com> Date: Mon, 8 Apr 2024 15:54:51 +0800 Subject: [PATCH 1/2] Can choose not to add dummy point in perodic case - This behavior is controlled by macro "INTP_PERIODIC_NO_DUMMY_POINT" - Prefix all macros with "INTP_", and rename some of them --- CMakeLists.txt | 7 +- src/include/BSpline.hpp | 22 ++-- src/include/BandMatrix.hpp | 20 ++-- src/include/DedicatedThreadPool.hpp | 18 ++-- src/include/Interpolation.hpp | 70 ++++++++----- src/include/InterpolationTemplate.hpp | 82 +++++++++------ src/include/Mesh.hpp | 8 ++ src/include/util.hpp | 16 ++- test/CMakeLists.txt | 14 ++- test/src/interpolation-speed-test.cpp | 36 ++++++- test/src/interpolation-template-test.cpp | 43 +++++++- test/src/interpolation-test.cpp | 124 +++++++++++++++++++---- 12 files changed, 339 insertions(+), 121 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 980323f..f6f0b64 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/src/include/BSpline.hpp b/src/include/BSpline.hpp index 663aab4..a653cd5 100644 --- a/src/include/BSpline.hpp +++ b/src/include/BSpline.hpp @@ -6,11 +6,10 @@ #include // fmod #include // ref #include // distance -#include // range_error #include // is_same, is_arithmatic #include -#ifdef _DEBUG +#ifdef INTP_DEBUG #include #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(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"; diff --git a/src/include/BandMatrix.hpp b/src/include/BandMatrix.hpp index b29f721..38f7e75 100644 --- a/src/include/BandMatrix.hpp +++ b/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 { } 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_); } diff --git a/src/include/DedicatedThreadPool.hpp b/src/include/DedicatedThreadPool.hpp index bfc8bef..62e0b6b 100644 --- a/src/include/DedicatedThreadPool.hpp +++ b/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 diff --git a/src/include/Interpolation.hpp b/src/include/Interpolation.hpp index 97a9786..ea8edfe 100644 --- a/src/include/Interpolation.hpp +++ b/src/include/Interpolation.hpp @@ -78,16 +78,21 @@ class InterpolationFunction { // TODO: Add integration DimArray&, std::pair 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(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 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(i) - @@ -120,23 +125,30 @@ class InterpolationFunction { // TODO: Add integration DimArray& input_coords, std::pair x_range) { uniform_[dim_ind] = false; + const auto periodic = periodicity(dim_ind); + const size_type n{static_cast( 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 +195,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 +461,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& @@ -469,7 +480,9 @@ class InterpolationFunction { // TODO: Add integration return spline_; } - size_type order() const { return order_; } + size_type order() const { + return order_; + } }; template @@ -485,10 +498,17 @@ class InterpolationFunction1D : public InterpolationFunction { : InterpolationFunction1D( std::make_pair(typename base::coord_type{}, static_cast( - (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 InterpolationFunction1D(std::pair x_range, diff --git a/src/include/InterpolationTemplate.hpp b/src/include/InterpolationTemplate.hpp index 0f7749f..30e0793 100644 --- a/src/include/InterpolationTemplate.hpp +++ b/src/include/InterpolationTemplate.hpp @@ -3,15 +3,18 @@ #include "BSpline.hpp" #include "BandLU.hpp" -#include "DedicatedThreadPool.hpp" #include "Mesh.hpp" #include "util.hpp" +#ifdef INTP_MULTITHREAD +#include "DedicatedThreadPool.hpp" +#endif + #if __cplusplus >= 201703L #include #endif -#ifdef _TRACE +#ifdef INTP_TRACE #include #endif @@ -231,6 +234,7 @@ class InterpolationFunctionTemplate { void build_solver_() { const auto& order = base_.order_; +#ifndef INTP_PERIODIC_NO_DUMMY_POINT // adjust dimension according to periodicity { DimArray dim_size_tmp; @@ -240,7 +244,7 @@ class InterpolationFunctionTemplate { } mesh_dimension_.resize(dim_size_tmp); } - +#endif DimArray base_spline_vals_per_dim; @@ -257,7 +261,7 @@ class InterpolationFunctionTemplate { } } -#ifdef _TRACE +#ifdef INTP_TRACE std::cout << "\n[TRACE] Coefficient Matrices\n"; #endif @@ -286,7 +290,7 @@ class InterpolationFunctionTemplate { mat_dim, band_width, band_width); #endif -#ifdef _TRACE +#ifdef INTP_TRACE std::cout << "\n[TRACE] Dimension " << d << '\n'; std::cout << "[TRACE] {0, 0} -> 1\n"; #endif @@ -384,7 +388,7 @@ class InterpolationFunctionTemplate { base_spline_vals_per_dim[d][j]; } #endif -#ifdef _TRACE +#ifdef INTP_TRACE std::cout << "[TRACE] {" << (periodic @@ -397,7 +401,7 @@ class InterpolationFunctionTemplate { } } -#ifdef _TRACE +#ifdef INTP_TRACE std::cout << "[TRACE] {" << mesh_dimension_.dim_size(d) - 1 << ", " << mesh_dimension_.dim_size(d) - 1 << "} -> 1\n"; #endif @@ -426,31 +430,47 @@ class InterpolationFunctionTemplate { ctrl_pt_type solve_for_control_points_( const Mesh& f_mesh) const { ctrl_pt_type weights{mesh_dimension_}; - ctrl_pt_type weights_tmp(1); // auxilary weight for swapping between - if CPP17_CONSTEXPR_ (dim > 1) { weights_tmp.resize(mesh_dimension_); } - - auto check_idx = - [&](typename Mesh::index_type& indices) { - bool keep_flag = true; - for (size_type d = 0; d < dim; ++d) { - if (base_.periodicity(d)) { - // Skip last point of periodic dimension - keep_flag = - keep_flag && indices[d] != weights.dim_size(d); - indices[d] = (indices[d] + weights.dim_size(d) + - base_.order_ / 2) % - weights.dim_size(d); - } - } - return keep_flag; - }; - - // Copy interpolating values into weights mesh as the initial state of - // the iterative control points solving algorithm +#ifdef INTP_PERIODIC_NO_DUMMY_POINT for (auto it = f_mesh.begin(); it != f_mesh.end(); ++it) { auto f_indices = f_mesh.iter_indices(it); - if (check_idx(f_indices)) { weights(f_indices) = *it; } + for (size_type d = 0; d < dim; ++d) { + if (base_.periodicity(d)) { + f_indices[d] = (f_indices[d] + weights.dim_size(d) + + base_.order_ / 2) % + weights.dim_size(d); + } + } + weights(f_indices) = *it; } +#else + { + auto check_idx = + [&](typename Mesh::index_type& indices) { + bool keep_flag = true; + for (size_type d = 0; d < dim; ++d) { + if (base_.periodicity(d)) { + // Skip last point of periodic dimension + keep_flag = + keep_flag && indices[d] != weights.dim_size(d); + indices[d] = (indices[d] + weights.dim_size(d) + + base_.order_ / 2) % + weights.dim_size(d); + } + } + return keep_flag; + }; + + // Copy interpolating values into weights mesh as the initial state + // of the iterative control points solving algorithm + for (auto it = f_mesh.begin(); it != f_mesh.end(); ++it) { + auto f_indices = f_mesh.iter_indices(it); + if (check_idx(f_indices)) { weights(f_indices) = *it; } + } + } +#endif + + ctrl_pt_type weights_tmp(1); // auxilary weight for swapping between + if CPP17_CONSTEXPR_ (dim > 1) { weights_tmp.resize(mesh_dimension_); } auto array_right_shift = [](DimArray arr) { DimArray new_arr{}; @@ -506,7 +526,7 @@ class InterpolationFunctionTemplate { } }; -#ifdef _MULTITHREAD +#ifdef INTP_MULTITHREAD // TODO: use a more robust task division strategy const size_type block_num = static_cast( std::sqrt(static_cast(hyperplane_size))); @@ -531,7 +551,7 @@ class InterpolationFunctionTemplate { for (auto&& f : res) { f.get(); } #else solve_and_rearrange_block(0, hyperplane_size); -#endif // _MULTITHREAD +#endif // INTP_MULTITHREAD } if CPP17_CONSTEXPR_ (dim % 2 == 0 || dim == 1) { diff --git a/src/include/Mesh.hpp b/src/include/Mesh.hpp index 2cab71c..c331a8f 100644 --- a/src/include/Mesh.hpp +++ b/src/include/Mesh.hpp @@ -260,6 +260,14 @@ class Mesh { const allocator_type& alloc = allocator_type()) : Mesh(std::make_pair(array.begin(), array.end()), alloc) {} + // convert constructor from mesh using different allocator + template + Mesh(const Mesh& mesh, + const allocator_type& alloc = allocator_type()) + : Mesh(mesh.dimension(), alloc) { + std::copy(mesh.begin(), mesh.end(), storage_.begin()); + } + // properties size_type size() const { return storage_.size(); } diff --git a/src/include/util.hpp b/src/include/util.hpp index 9de4999..c1abbf8 100644 --- a/src/include/util.hpp +++ b/src/include/util.hpp @@ -67,7 +67,7 @@ void dispatch_indexed(Func&& func, Args&&... args) { std::forward(args)...); } -#ifdef STACK_ALLOCATOR +#ifdef INTP_STACK_ALLOCATOR /** * @brief A simple stack allocator, with fixed size and a LIFO allocation @@ -244,11 +244,17 @@ struct lazy_conditional { using type = FalseTemplate; }; -#ifdef _DEBUG -#define CUSTOM_ASSERT(assertion, msg) \ - if (!(assertion)) { throw std::runtime_error(msg); } +#ifdef INTP_DEBUG +#define INTP_ENABLE_ASSERTION +#endif + +#ifdef INTP_ENABLE_ASSERTION +#define INTP_ASSERT(assertion, msg) \ + do { \ + if (!(assertion)) { throw std::runtime_error(msg); } \ + } while (0) #else -#define CUSTOM_ASSERT(assertion, msg) +#define INTP_ASSERT(assertion, msg) #endif /** diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 47f7524..9c52bdd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,6 +3,12 @@ cmake_minimum_required(VERSION 3.12.0) set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_EXTENSIONS OFF) +# Set default build type to debug. +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Debug) +endif() +message(STATUS "Build type = " ${CMAKE_BUILD_TYPE}) + # Add options according to compiler. if(MSVC) set(CMAKE_CXX_FLAGS "/W4 /EHsc") @@ -18,9 +24,11 @@ else() endif() endif() -add_compile_definitions($<$:_TRACE>) -add_compile_definitions($<$,$>:_DEBUG>) -add_compile_definitions($<$,$>:_MULTITHREAD>) +add_compile_definitions($<$:INTP_TRACE> + $<$,$>:INTP_DEBUG> + $<$,$>:INTP_MULTITHREAD>) + +add_compile_definitions(INTP_PERIODIC_NO_DUMMY_POINT) # Specify tests list(APPEND tests "util-test" "mesh-test" "band-matrix-and-solver-test" "bspline-test" "interpolation-test" "interpolation-speed-test" "interpolation-template-test") diff --git a/test/src/interpolation-speed-test.cpp b/test/src/interpolation-speed-test.cpp index 69dfb9c..81e3150 100644 --- a/test/src/interpolation-speed-test.cpp +++ b/test/src/interpolation-speed-test.cpp @@ -62,12 +62,20 @@ int main() { constexpr size_t len_1d = len * len * len; constexpr double dx = 2 * M_PI / (len_1d); std::vector vec_1d{}; + +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + vec_1d.reserve(len_1d); + for (size_t i = 0; i < len_1d; ++i) { + vec_1d.emplace_back( + std::sin(13 * (static_cast(i) * dx - M_PI))); + } +#else vec_1d.reserve(len_1d + 1); for (size_t i = 0; i <= len_1d; ++i) { vec_1d.emplace_back( std::sin(13 * (static_cast(i) * dx - M_PI))); } - +#endif std::vector vals_1d; vals_1d.reserve(coord_1d.size()); for (auto& x : coord_1d) { vals_1d.push_back(std::sin(13 * x)); } @@ -121,6 +129,16 @@ int main() { const size_t len_2d = static_cast(std::pow(len, 1.5)); const double dt = 2 * M_PI / static_cast(len_2d); +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh trig_mesh_2d(len_2d); + for (size_t i = 0; i < len_2d; ++i) { + for (size_t j = 0; j < len_2d; ++j) { + trig_mesh_2d(i, j) = + std::sin(static_cast(i) * dt - M_PI) * + std::cos(static_cast(j) * dt - M_PI); + } + } +#else Mesh trig_mesh_2d(len_2d + 1); for (size_t i = 0; i <= len_2d; ++i) { for (size_t j = 0; j <= len_2d; ++j) { @@ -129,6 +147,7 @@ int main() { std::cos(static_cast(j) * dt - M_PI); } } +#endif std::vector vals_2d; vals_2d.reserve(coord_2d.size()); @@ -190,6 +209,20 @@ int main() { constexpr double dt_3d = 2 * M_PI / len; constexpr double dt_3d_aperiodic = 1. / (len - 1); +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh mesh_3d{len, len, len}; + for (size_t i = 0; i < len; ++i) { + for (size_t j = 0; j < len; ++j) { + for (size_t k = 0; k < len; ++k) { + mesh_3d(i, j, k) = + std::sin(static_cast(i) * dt_3d - M_PI) * + std::cos(static_cast(j) * dt_3d - M_PI) * + std::exp(-std::pow( + static_cast(k) * dt_3d_aperiodic - .5, 2)); + } + } + } +#else Mesh mesh_3d{len + 1, len + 1, len}; for (size_t i = 0; i <= len; ++i) { for (size_t j = 0; j <= len; ++j) { @@ -202,6 +235,7 @@ int main() { } } } +#endif std::vector vals_3d; vals_3d.reserve(coord_3d.size()); diff --git a/test/src/interpolation-template-test.cpp b/test/src/interpolation-template-test.cpp index dd45e2c..37e217f 100644 --- a/test/src/interpolation-template-test.cpp +++ b/test/src/interpolation-template-test.cpp @@ -48,11 +48,19 @@ int main() { constexpr size_t len = 1024; constexpr double dx = 2 * M_PI / (len * len); std::vector trig_vec{}; +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + trig_vec.reserve(len * len); + for (size_t i = 0; i < len * len; ++i) { + trig_vec.emplace_back( + std::sin(7 * (static_cast(i) * dx - M_PI))); + } +#else trig_vec.reserve(len * len + 1); for (size_t i = 0; i < len * len + 1; ++i) { trig_vec.emplace_back( std::sin(7 * (static_cast(i) * dx - M_PI))); } +#endif std::vector vals_1d; vals_1d.reserve(coord_1d.size()); @@ -106,6 +114,20 @@ int main() { constexpr double dt = 2 * M_PI / len; +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh trig_mesh_2d_1(len); + Mesh trig_mesh_2d_2(len); + for (size_t i = 0; i < len; ++i) { + for (size_t j = 0; j < len; ++j) { + trig_mesh_2d_1(i, j) = + std::sin(static_cast(i) * dt - M_PI) * + std::cos(static_cast(j) * dt - M_PI); + trig_mesh_2d_2(i, j) = + std::cos(2 * (static_cast(i) * dt - M_PI)) * + std::sin(2 * (static_cast(j) * dt - M_PI)); + } + } +#else Mesh trig_mesh_2d_1(len + 1); Mesh trig_mesh_2d_2(len + 1); for (size_t i = 0; i <= len; ++i) { @@ -118,7 +140,7 @@ int main() { std::sin(2 * (static_cast(j) * dt - M_PI)); } } - +#endif std::vector vals_2d_1; std::vector vals_2d_2; vals_2d_1.reserve(coord_2d.size()); @@ -187,6 +209,24 @@ int main() { constexpr double dt_3d = 2 * M_PI / lt; constexpr double dz = 1. / (lz - 1); +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh mesh_3d_1(lt, lt, lz); + Mesh mesh_3d_2(lt, lt, lz); + for (size_t i = 0; i < lt; ++i) { + for (size_t j = 0; j < lt; ++j) { + for (size_t k = 0; k < lz; ++k) { + mesh_3d_1(i, j, k) = + std::sin(static_cast(i) * dt_3d - M_PI) * + std::cos(static_cast(j) * dt_3d - M_PI) * + std::exp(-std::pow(static_cast(k) * dz - .5, 2)); + mesh_3d_2(i, j, k) = + std::cos(2 * (static_cast(i) * dt_3d - M_PI)) * + std::sin(2 * (static_cast(j) * dt_3d - M_PI)) * + std::exp(-std::pow(static_cast(k) * dz - .5, 3)); + } + } + } +#else Mesh mesh_3d_1(lt + 1, lt + 1, lz); Mesh mesh_3d_2(lt + 1, lt + 1, lz); for (size_t i = 0; i <= lt; ++i) { @@ -203,6 +243,7 @@ int main() { } } } +#endif std::vector vals_3d_1; std::vector vals_3d_2; diff --git a/test/src/interpolation-test.cpp b/test/src/interpolation-test.cpp index 97f692f..868eacb 100644 --- a/test/src/interpolation-test.cpp +++ b/test/src/interpolation-test.cpp @@ -61,14 +61,18 @@ int main() { util::get_range(vals_1d_linear)); assertion(d < tol); std::cout << "\n1D test (linear) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; d = rel_err(interp1, util::get_range(coords_1d_half), util::get_range(vals_1d)); assertion(d < tol); std::cout << "\n1D test (cubic) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; assertion(std::abs(interp1(-.5) - (-6.3167718907512755)) < tol, @@ -128,17 +132,19 @@ int main() { 0.2755974099701995}}; assertion(interp2.uniform(0) && interp2.uniform(1), - "Uniform properties check failed."); + "Uniform properties check \033[1;91mfailed\033[0m."); d = rel_err(interp2, util::get_range(coords_2d), util::get_range(vals_2d)); assertion(d < tol); std::cout << "\n2D test " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; try { interp2.at(-1, 1); - assertion(false, "Out of boundary check failed.\n"); + assertion(false, "Out of boundary check \033[1;91mfailed\033[0m.\n"); } catch (const std::domain_error&) { std::cout << "Out of boundary check succeed.\n"; } @@ -277,7 +283,9 @@ int main() { d = rel_err(interp3, util::get_range(coords_3d), util::get_range(vals_3d)); assertion(d < tol); std::cout << "\n3D test " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; assertion(!interp3.periodicity(0) && !interp3.periodicity(1) && @@ -287,6 +295,14 @@ int main() { std::cout << "\n1D Interpolation with Periodic Boundary Test:\n"; +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + const auto f_periodic_range = std::make_pair(f.begin(), std::prev(f.end())); + InterpolationFunction1D<> interp1_periodic_linear(f_periodic_range, 1, + true); + + InterpolationFunction1D<> interp1_periodic(f_periodic_range, 4, true); + +#else // Since the InterpolationFunction ignores the last term along periodic // dimension, thus we can use the origin non-periodic data to interpolate a // periodic spline function @@ -295,6 +311,7 @@ int main() { true); InterpolationFunction1D<> interp1_periodic(util::get_range(f), 4, true); +#endif // linear interpolation results std::vector vals_1d_periodic_linear; @@ -316,14 +333,18 @@ int main() { util::get_range(vals_1d_periodic_linear)); assertion(d < tol); std::cout << "\n1D test (linear) with periodic boundary " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; d = rel_err(interp1_periodic, util::get_range(coords_1d), util::get_range(vals_1d_periodic)); assertion(d < tol); std::cout << "\n1D test (quartic) with periodic boundary " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; assertion(std::abs(interp1_periodic(*coords_1d.begin() - (f.size() - 1)) - @@ -336,10 +357,24 @@ int main() { // 2D interpolation test with one dimension being periodic boundary std::cout << "\n2D Interpolation with one Periodic Boundary Test:\n"; +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh f2d_periodic_y{f2.size(), f2[0].size() - 1}; + for (size_t i = 0; i < f2d_periodic_y.dim_size(0); ++i) { + for (size_t j = 0; j < f2d_periodic_y.dim_size(1); ++j) { + f2d_periodic_y(i, j) = f2[i][j]; + } + } + InterpolationFunction interp2_periodic( + 3, {false, true}, f2d_periodic_y, + make_pair(0., static_cast(f2d.dim_size(0)) - 1.), + make_pair(0., static_cast(f2d.dim_size(1)) - 1.)); + +#else InterpolationFunction interp2_periodic( 3, {false, true}, f2d, make_pair(0., static_cast(f2d.dim_size(0)) - 1.), make_pair(0., static_cast(f2d.dim_size(1)) - 1.)); +#endif auto vals_2d_periodic = {-0.5456439415470818, 0.7261218483070795, 0.21577722210958022, -1.6499933881987376, @@ -351,7 +386,9 @@ int main() { util::get_range(vals_2d_periodic)); assertion(d < tol); std::cout << "\n2D test with periodic boundary " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; #ifdef _DEBUG @@ -375,7 +412,9 @@ int main() { util::get_range(coords_1d_half), util::get_range(vals_1d_derivative_1)); assertion(d < tol); std::cout << "\n1D test of derivative " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; auto vals_1d_derivative_periodic = { @@ -392,7 +431,9 @@ int main() { util::get_range(vals_1d_derivative_periodic)); assertion(d < tol); std::cout << "\n1D test of derivative in periodic case " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; // 2D interpolation derivative test @@ -412,7 +453,9 @@ int main() { util::get_range(coords_2d), util::get_range(vals_2d_derivative_x2_y1)); assertion(d < tol); std::cout << "\n2D test of derivative (2,1) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; // 3D interpolation derivative test @@ -433,7 +476,9 @@ int main() { util::get_range(vals_3d_derivative_x1_y0_z3)); assertion(d < tol); std::cout << "\n3D test of derivative (1,0,3) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; // 1D non-uniform interpolation test @@ -483,25 +528,37 @@ int main() { util::get_range(vals_1d_nonuniform_linear)); assertion(d < tol); std::cout << "\n1D nonuniform test (linear) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; d = rel_err(interp1_nonuniform, util::get_range(coords_1d), util::get_range(vals_1d_nonuniform)); assertion(d < tol); std::cout << "\n1D nonuniform test (cubic) " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; std::cout << "\n1D nonuniform Interpolation Test with Periodic Boundary:\n"; // 1D non-uniform periodic interpolation test +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + InterpolationFunction1D<> interp1_nonuniform_periodic_linear( + util::get_range(input_coords_1d), f_periodic_range, 1, true); + + InterpolationFunction1D<> interp1_nonuniform_periodic( + util::get_range(input_coords_1d), f_periodic_range, 4, true); +#else InterpolationFunction1D<> interp1_nonuniform_periodic_linear( util::get_range(input_coords_1d), util::get_range(f), 1, true); InterpolationFunction1D<> interp1_nonuniform_periodic( util::get_range(input_coords_1d), util::get_range(f), 4, true); +#endif // linear interpolation results std::vector vals_1d_nonuniform_periodic_linear; @@ -525,12 +582,19 @@ int main() { d = rel_err(interp1_nonuniform_periodic_linear, util::get_range(coords_1d), util::get_range(vals_1d_nonuniform_periodic_linear)); assertion(d < tol); + std::cout << "\n1D nonuniform test (linear) with periodic boundary " + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; + std::cout << "Relative Error = " << d << '\n'; d = rel_err(interp1_nonuniform_periodic, util::get_range(coords_1d), util::get_range(vals_1d_nonuniform_periodic)); assertion(d < tol); - std::cout << "\n1D nonuniform test with periodic boundary " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + std::cout << "\n1D nonuniform test (quartic) with periodic boundary " + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; // 2D non-uniform with one dimension being periodic boundary @@ -539,10 +603,24 @@ int main() { auto nonuniform_coord_for_2d = {0., 1.329905262345947, 2.200286645200226, 3.1202827237815516, 4.}; + +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + Mesh f2d_periodic_x{f2.size() - 1, f2[0].size()}; + for (size_t i = 0; i < f2d_periodic_x.dim_size(0); ++i) { + for (size_t j = 0; j < f2d_periodic_x.dim_size(1); ++j) { + f2d_periodic_x(i, j) = f2[i][j]; + } + } + InterpolationFunction interp2_X_periodic_Y_nonuniform( + 3, {true, false}, f2d_periodic_x, + std::make_pair(0., static_cast(f2d.dim_size(0)) - 1.), + util::get_range(nonuniform_coord_for_2d)); +#else InterpolationFunction interp2_X_periodic_Y_nonuniform( 3, {true, false}, f2d, std::make_pair(0., static_cast(f2d.dim_size(0)) - 1.), util::get_range(nonuniform_coord_for_2d)); +#endif auto vals_2d_X_periodic_Y_nonuniform = { -0.7160424002258807, 0.6702846215275077, 0.04933393138376427, @@ -554,13 +632,20 @@ int main() { util::get_range(vals_2d_X_periodic_Y_nonuniform)); assertion(d < tol); std::cout << "\n2D test with x-periodic and y-nonuniform " - << (assertion.last_status() == 0 ? "succeed" : "failed") << '\n'; + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") + << '\n'; std::cout << "Relative Error = " << d << '\n'; // interpolation on vector + std::vector> pts; constexpr std::size_t pts_n = 31; +#ifdef INTP_PERIODIC_NO_DUMMY_POINT + for (std::size_t i = 0; i < pts_n; ++i) { +#else for (std::size_t i = 0; i <= pts_n; ++i) { +#endif double theta = 2. * M_PI * static_cast(i) / pts_n; pts.emplace_back(std::cos(theta), std::sin(theta)); } @@ -579,7 +664,8 @@ int main() { err /= sample_n; assertion(err < 1.e-5f); std::cout << "Interpolation of 2d points on a circle " - << (assertion.last_status() == 0 ? "succeed" : "failed") + << (assertion.last_status() == 0 ? "succeed" + : "\033[1;91mfailed\033[0m") << '\n'; std::cout << "Relative Error = " << err << '\n'; } From fa1fcaaf4170997998f96d6e81b18ccca62adfb6 Mon Sep 17 00:00:00 2001 From: qzhong <12ff54e@gmail.com> Date: Mon, 8 Apr 2024 19:43:43 +0800 Subject: [PATCH 2/2] Add description of macros in README --- CMakeLists.txt | 2 +- README.md | 14 +++++++++++++- src/include/Interpolation.hpp | 4 ++++ 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f6f0b64..aa45bf8 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/README.md b/README.md index e9738b7..614884c 100644 --- a/README.md +++ b/README.md @@ -42,6 +42,18 @@ InterpolationFunction 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 diff --git a/src/include/Interpolation.hpp b/src/include/Interpolation.hpp index ea8edfe..41a9499 100644 --- a/src/include/Interpolation.hpp +++ b/src/include/Interpolation.hpp @@ -4,6 +4,10 @@ #include // ceil #include +#ifdef INTP_TRACE +#define INTP_DEBUG +#endif + #include "InterpolationTemplate.hpp" namespace intp {