add predefined triangulation, -Ofast compilation

.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"
+    }
+}

imops/cpp/src/interp2d/interpolator.h

@@ -3,46 +3,48 @@
 class Interpolator {
 
 public:
-    using pyarr_float = py::array_t<float, py::array::c_style | py::array::forcecast>;
+    using pyarr_double = py::array_t<double, py::array::c_style | py::array::forcecast>;
     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<Triangulator::pyarr_size_t> 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<float> int_values(n, fill_value);
+        std::vector<double> 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
         for (size_t i = 0; i < n; ++i) {
             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<double>(coords_info[3]);
-                double lambda_1 = static_cast<double>(coords_info[0]) * one_over_2area;
-                double lambda_2 = static_cast<double>(coords_info[1]) * one_over_2area;
-                double lambda_3 = static_cast<double>(coords_info[2]) * one_over_2area;
-                double f1 = values.at(triangulation.triangles.at(t));  // implicit float -> double
+                double lambda_1 = static_cast<double>(coords_info[0]);
+                double lambda_2 = static_cast<double>(coords_info[1]);
+                double lambda_3 = static_cast<double>(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<float>(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()};
     }
 };

imops/cpp/src/interp2d/triangulator.h

@@ -4,6 +4,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 #include <pybind11/embed.h>
+#include <pybind11/stl.h>
 #include "delaunator/delaunator-header-only.hpp"
 #include "utils.h"
 
@@ -21,31 +22,57 @@ class Triangulator {
     std::vector<std::vector<size_t>> point2tri;
     std::unordered_map<uint64_t, std::array<size_t, 2>> edge2tri;
 
-    Triangulator(const pyarr_size_t& pypoints, int n_jobs = 1) {
+    Triangulator(const pyarr_size_t& pypoints, int n_jobs, std::optional<pyarr_size_t> 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> 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<double>(pypoints.at(i, 0));
-            double_points.at(2 * i + 1) = static_cast<double>(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<size_t>(delaunated.coords.at(j));
-            points.at(j + 1) = static_cast<size_t>(delaunated.coords.at(j + 1));
+
+        else {
+            std::vector<double> double_points(2 * n);
+            #pragma parallel for
+            for (size_t i = 0; i < n; ++i) {
+                double_points.at(2 * i) = static_cast<double>(pypoints.at(i, 0));
+                double_points.at(2 * i + 1) = static_cast<double>(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<size_t>(delaunated.coords.at(j));
+                points.at(j + 1) = static_cast<size_t>(delaunated.coords.at(j + 1));
+            }
         }
+
         point2tri.resize(n);
         for (size_t i = 0; i < triangles.size() / 3; ++i) {
             size_t t = 3 * i;

imops/cpp/src/main.cpp

@@ -5,7 +5,7 @@
 
 PYBIND11_MODULE(cpp_modules, m) {
     py::class_<Interpolator>(m, "Linear2DInterpolatorCpp", "Interpolator class")
-        .def(py::init<const Triangulator::pyarr_size_t&, int>())
+        .def(py::init<const Triangulator::pyarr_size_t&, int, std::optional<Triangulator::pyarr_size_t>>())
         .def("__call__", &Interpolator::operator());
 }
 

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