From 3567f1ced4dad239860f377f5f78610b4c8163cd Mon Sep 17 00:00:00 2001 From: Alexey Belkov Date: Sat, 23 Dec 2023 02:13:12 +0300 Subject: [PATCH] add predefined triangulation, -Ofast compilation --- .vscode/settings.json | 54 +++++++++++++++++++++++++ imops/cpp/src/interp2d/interpolator.h | 34 ++++++++-------- imops/cpp/src/interp2d/triangulator.h | 57 ++++++++++++++++++++------- imops/cpp/src/main.cpp | 2 +- imops/interp2d.py | 4 +- 5 files changed, 117 insertions(+), 34 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 00000000..33755e1c --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,54 @@ +{ + "files.associations": { + "*.config": "python", + "array": "cpp", + "atomic": "cpp", + "bit": "cpp", + "*.tcc": "cpp", + "cctype": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "compare": "cpp", + "complex": "cpp", + "concepts": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdint": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "deque": "cpp", + "string": "cpp", + "unordered_map": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "functional": "cpp", + "iterator": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "fstream": "cpp", + "initializer_list": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "limits": "cpp", + "new": "cpp", + "numbers": "cpp", + "ostream": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "streambuf": "cpp", + "cinttypes": "cpp", + "typeinfo": "cpp" + } +} \ No newline at end of file diff --git a/imops/cpp/src/interp2d/interpolator.h b/imops/cpp/src/interp2d/interpolator.h index 68ebea17..b816134b 100644 --- a/imops/cpp/src/interp2d/interpolator.h +++ b/imops/cpp/src/interp2d/interpolator.h @@ -3,25 +3,27 @@ class Interpolator { public: - using pyarr_float = py::array_t; + using pyarr_double = py::array_t; const Triangulator triangulation; - Interpolator(const Triangulator::pyarr_size_t& points, int n_jobs) - : triangulation(points, n_jobs) { + Interpolator(const Triangulator::pyarr_size_t& pypoints, int n_jobs, std::optional pytriangles) + : triangulation(pypoints, n_jobs, pytriangles) { + // if (pytriangles.has_value()) + // std::cout << pytriangles -> shape(0) << ' ' << pytriangles -> shape(1) << '\n'; } - pyarr_float operator()(const Triangulator::pyarr_size_t& int_points, const pyarr_float& values, + pyarr_double operator()(const Triangulator::pyarr_size_t& int_points, const pyarr_double& values, const Triangulator::pyarr_size_t& neighbors, double fill_value = 0.0) { - if (triangulation.points.size() / 2 != values.shape()[0]) { + if (triangulation.points.size() / 2 != values.shape(0)) { throw std::invalid_argument("Length mismatch between known points and their values"); } - if (neighbors.shape()[0] != int_points.shape()[0]) { + if (neighbors.shape(0) != int_points.shape(0)) { throw std::invalid_argument("Length mismatch between int_points and their neighbors"); } size_t n = int_points.shape()[0]; - std::vector int_values(n, fill_value); + std::vector int_values(n, fill_value); - omp_set_dynamic(0); // Explicitly disable dynamic teams + omp_set_dynamic(0); // Explicitly disable dynamic teams omp_set_num_threads(triangulation.n_jobs_); #pragma parallel for @@ -29,20 +31,20 @@ class Interpolator { size_t x = int_points.at(i, 0); size_t y = int_points.at(i, 1); auto location_info = triangulation.locate_point(x, y, neighbors.at(i)); - if (location_info != std::nullopt) { + if (location_info.has_value()) { size_t t = location_info->first; auto& coords_info = location_info->second; double one_over_2area = 1.0 / static_cast(coords_info[3]); - double lambda_1 = static_cast(coords_info[0]) * one_over_2area; - double lambda_2 = static_cast(coords_info[1]) * one_over_2area; - double lambda_3 = static_cast(coords_info[2]) * one_over_2area; - double f1 = values.at(triangulation.triangles.at(t)); // implicit float -> double + double lambda_1 = static_cast(coords_info[0]); + double lambda_2 = static_cast(coords_info[1]); + double lambda_3 = static_cast(coords_info[2]); + double f1 = values.at(triangulation.triangles.at(t)); double f2 = values.at(triangulation.triangles.at(t + 1)); double f3 = values.at(triangulation.triangles.at(t + 2)); - int_values[i] = static_cast(f1 * lambda_1 + f2 * lambda_2 + f3 * lambda_3); + int_values[i] = (f1 * lambda_1 + f2 * lambda_2 + f3 * lambda_3) * one_over_2area; } } - - return {{n}, {sizeof(float)}, int_values.data()}; + + return {{n}, {sizeof(double)}, int_values.data()}; } }; diff --git a/imops/cpp/src/interp2d/triangulator.h b/imops/cpp/src/interp2d/triangulator.h index ff04210e..1644e553 100644 --- a/imops/cpp/src/interp2d/triangulator.h +++ b/imops/cpp/src/interp2d/triangulator.h @@ -4,6 +4,7 @@ #include #include #include +#include #include "delaunator/delaunator-header-only.hpp" #include "utils.h" @@ -21,31 +22,57 @@ class Triangulator { std::vector> point2tri; std::unordered_map> edge2tri; - Triangulator(const pyarr_size_t& pypoints, int n_jobs = 1) { + Triangulator(const pyarr_size_t& pypoints, int n_jobs, std::optional pytriangles) { if (n_jobs < 0 and n_jobs != -1) { throw std::invalid_argument( "Invalid number of workers, have to be -1 or positive integer"); } n_jobs_ = n_jobs == -1 ? omp_get_num_procs() : n_jobs; - size_t n = pypoints.shape()[0]; - std::vector double_points(2 * n); + size_t n = pypoints.shape(0); + points.resize(2 * n); + omp_set_dynamic(0); // Explicitly disable dynamic teams omp_set_num_threads(n_jobs_); - #pragma parallel for - for (size_t i = 0; i < n; ++i) { - double_points.at(2 * i) = static_cast(pypoints.at(i, 0)); - double_points.at(2 * i + 1) = static_cast(pypoints.at(i, 1)); + + if (pytriangles.has_value()) { + #pragma parallel for + for (size_t i = 0; i < n; ++i) { + size_t j = 2 * i; + points.at(j) = pypoints.at(i, 0); + points.at(j + 1) = pypoints.at(i, 1); + } + + size_t m = pytriangles -> shape(0); + triangles.resize(3 * m); + + #pragma parallel for + for (size_t i = 0; i < m; ++i) { + size_t j = 3 * i; + triangles.at(j) = pytriangles -> at(i, 0); + triangles.at(j + 1) = pytriangles -> at(i, 1); + triangles.at(j + 2) = pytriangles -> at(i, 2); + } } - delaunator::Delaunator delaunated(double_points); - triangles = std::move(delaunated.triangles); - - #pragma parallel for - for (size_t i = 0; i < n; ++i) { - size_t j = 2 * i; - points.at(j) = static_cast(delaunated.coords.at(j)); - points.at(j + 1) = static_cast(delaunated.coords.at(j + 1)); + + else { + std::vector double_points(2 * n); + #pragma parallel for + for (size_t i = 0; i < n; ++i) { + double_points.at(2 * i) = static_cast(pypoints.at(i, 0)); + double_points.at(2 * i + 1) = static_cast(pypoints.at(i, 1)); + } + delaunator::Delaunator delaunated(double_points); + triangles = std::move(delaunated.triangles); + + #pragma parallel for + for (size_t i = 0; i < n; ++i) { + size_t j = 2 * i; + points.at(j) = static_cast(delaunated.coords.at(j)); + points.at(j + 1) = static_cast(delaunated.coords.at(j + 1)); + } } + point2tri.resize(n); for (size_t i = 0; i < triangles.size() / 3; ++i) { size_t t = 3 * i; diff --git a/imops/cpp/src/main.cpp b/imops/cpp/src/main.cpp index 585b9e4d..25817d83 100644 --- a/imops/cpp/src/main.cpp +++ b/imops/cpp/src/main.cpp @@ -5,7 +5,7 @@ PYBIND11_MODULE(cpp_modules, m) { py::class_(m, "Linear2DInterpolatorCpp", "Interpolator class") - .def(py::init()) + .def(py::init>()) .def("__call__", &Interpolator::operator()); } diff --git a/imops/interp2d.py b/imops/interp2d.py index 7717f310..8dc3c31b 100644 --- a/imops/interp2d.py +++ b/imops/interp2d.py @@ -22,8 +22,8 @@ class Linear2DInterpolator(Linear2DInterpolatorCpp): interp_values = interpolator(interp_points, values + 1.0, fill_value=0.0) """ - def __init__(self, points: np.ndarray, values: np.ndarray, n_jobs: int = 1, **kwargs): - super().__init__(points, n_jobs) + def __init__(self, points: np.ndarray, values: np.ndarray, n_jobs: int = 1, triangles: Optional[np.ndarray] = None, **kwargs): + super().__init__(points, n_jobs, triangles) self.kdtree = KDTree(data=points, **kwargs) self.values = values self.n_jobs = n_jobs