add test draft and test data

neuro-ml · Dec 24, 2023 · a47c6d3 · a47c6d3
1 parent 3567f1c
commit a47c6d3
Show file tree

Hide file tree

Showing 16 changed files with 72 additions and 27 deletions.
diff --git a/imops/cpp/src/interp2d/interpolator.h b/imops/cpp/src/interp2d/interpolator.h
@@ -5,14 +5,14 @@ class Interpolator {
 public:
     using pyarr_double = py::array_t<double, py::array::c_style | py::array::forcecast>;
     const Triangulator triangulation;
-    Interpolator(const Triangulator::pyarr_size_t& pypoints, int n_jobs, std::optional<Triangulator::pyarr_size_t> pytriangles)
+    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_double operator()(const Triangulator::pyarr_size_t& int_points, const pyarr_double& values,
-                           const Triangulator::pyarr_size_t& neighbors, double fill_value = 0.0) {
+    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)) {
             throw std::invalid_argument("Length mismatch between known points and their values");
         }
@@ -21,7 +21,7 @@ class Interpolator {
         }
 
         size_t n = int_points.shape()[0];
-        std::vector<double> int_values(n, fill_value);
+        std::vector<double> int_values(n);
 
         omp_set_dynamic(0);  // Explicitly disable dynamic teams
         omp_set_num_threads(triangulation.n_jobs_);
@@ -42,6 +42,8 @@ class Interpolator {
                 double f2 = values.at(triangulation.triangles.at(t + 1));
                 double f3 = values.at(triangulation.triangles.at(t + 2));
                 int_values[i] = (f1 * lambda_1 + f2 * lambda_2 + f3 * lambda_3) * one_over_2area;
+            } else {
+                int_values[i] = fill_value;
             }
         }
 

diff --git a/imops/cpp/src/interp2d/triangulator.h b/imops/cpp/src/interp2d/triangulator.h
@@ -22,7 +22,8 @@ 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, std::optional<pyarr_size_t> pytriangles) {
+    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");
@@ -42,16 +43,16 @@ class Triangulator {
                 points.at(j) = pypoints.at(i, 0);
                 points.at(j + 1) = pypoints.at(i, 1);
             }
-            
-            size_t m = pytriangles -> shape(0);
+
+            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);
+                triangles.at(j) = pytriangles->at(i, 0);
+                triangles.at(j + 1) = pytriangles->at(i, 1);
+                triangles.at(j + 2) = pytriangles->at(i, 2);
             }
         }
 

diff --git a/imops/interp2d.py b/imops/interp2d.py
@@ -7,20 +7,20 @@
 
 
 class Linear2DInterpolator(Linear2DInterpolatorCpp):
-    """
-    2D delaunay triangulation and parallel linear interpolation
-    Example:
-    n, m = 1024, 2
-    points = np.random.randint(low=0, high=1024, size=(n, m))
-    points = np.unique(points, axis=0)
-    x_points = points[: n // 2]
-    values = np.random.uniform(low=0.0, high=1.0, size=(len(x_points),))
-    interp_points = points[n // 2:]
-    n_jobs = -1 # will be equal to num of CPU cores
-    interpolator = Linear2DInterpolator(x_points, values, n_jobs)
-    # Also you can pass values to __call__ and rewrite the ones that were passed to __init__
-    interp_values = interpolator(interp_points, values + 1.0, fill_value=0.0)
-    """
+    # """
+    # #### 2D delaunay triangulation and parallel linear interpolation     
+    # Example:      
+    # n, m = 1024, 2        
+    # points = np.random.randint(low=0, high=1024, size=(n, m))        
+    # points = np.unique(points, axis=0)        
+    # x_points = points[: n // 2]        
+    # values = np.random.uniform(low=0.0, high=1.0, size=(len(x_points),))        
+    # interp_points = points[n // 2:]        
+    # n_jobs = -1 // will be equal to num of CPU cores        
+    # interpolator = Linear2DInterpolator(x_points, values, n_jobs)        
+    # // Also you can pass values to __call__ and rewrite the ones that were passed to __init__        
+    # 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,  triangles: Optional[np.ndarray] = None, **kwargs):
         super().__init__(points, n_jobs, triangles)

diff --git a/tests/test_data/arr_0.npy b/tests/test_data/arr_0.npy
diff --git a/tests/test_data/arr_1.npy b/tests/test_data/arr_1.npy
diff --git a/tests/test_data/arr_2.npy b/tests/test_data/arr_2.npy
diff --git a/tests/test_data/arr_3.npy b/tests/test_data/arr_3.npy
diff --git a/tests/test_data/arr_4.npy b/tests/test_data/arr_4.npy
diff --git a/tests/test_data/arr_5.npy b/tests/test_data/arr_5.npy
diff --git a/tests/test_data/arr_6.npy b/tests/test_data/arr_6.npy
diff --git a/tests/test_data/val_2.npy b/tests/test_data/val_2.npy
diff --git a/tests/test_data/val_3.npy b/tests/test_data/val_3.npy
diff --git a/tests/test_data/val_4.npy b/tests/test_data/val_4.npy
diff --git a/tests/test_data/val_5.npy b/tests/test_data/val_5.npy
diff --git a/tests/test_data/val_6.npy b/tests/test_data/val_6.npy
diff --git a/tests/test_interp2d.py b/tests/test_interp2d.py
@@ -0,0 +1,42 @@
+import numpy as np
+import scipy
+from scipy.interpolate import griddata
+
+from imops.interp2d import Linear2DInterpolator
+from deli import load_numpy
+
+np.random.seed(1337)
+
+
+def test_test_data():
+    for i in range(2):
+        x = load_numpy(f'test_data/arr_{i}.npy')
+        int_points = np.transpose((np.isnan(x)).nonzero())
+        x_points: np.ndarray = np.transpose((~np.isnan(x)).nonzero())
+        x_values = x[~np.isnan(x)]
+        imops_values = Linear2DInterpolator(x_points, x_values, n_jobs=1)(int_points, fill_value=0.0)
+        triangles = scipy.spatial.Delaunay(x_points).simplices
+        delaunay_values = Linear2DInterpolator(x_points, x_values, n_jobs=1, triangles=triangles)(int_points, fill_value=0.0)
+        scipy_values = griddata(x_points, x_values, int_points, method='linear', fill_value=0.0)
+
+        delta_ds = np.abs(delaunay_values - scipy_values)
+        # delta_di = np.abs(delaunay_values - imops_values)
+        delta_si = np.abs(scipy_values - imops_values)
+
+        assert delta_ds.max() <= 1e-8 and delta_si.max() <= 5, 'Failed with big cases, arr_0 or arr_1'
+
+    for i in range(2, 7):
+        x = load_numpy(f'test_data/arr_{i}.npy')
+        x_values = load_numpy(f'test_data/val_{i}.npy')
+        x_points = np.transpose(x.nonzero())
+        int_points = np.transpose((~x).nonzero())
+        imops_values = Linear2DInterpolator(x_points, x_values, n_jobs=1)(int_points, fill_value=0.0)
+        triangles = scipy.spatial.Delaunay(x_points).simplices
+        delaunay_values = Linear2DInterpolator(x_points, x_values, n_jobs=1, triangles=triangles)(int_points, fill_value=0.0)
+        scipy_values = griddata(x_points, x_values, int_points, method='linear', fill_value=0.0)
+
+        delta_ds = np.abs(delaunay_values - scipy_values)
+        delta_di = np.abs(delaunay_values - imops_values)
+        delta_si = np.abs(scipy_values - imops_values)
+
+        assert delta_di.max() <= 1.5 and delta_ds.max() <= 1e-10 and delta_si.max() <= 1.5, 'Failed with small cases, arr_2 ...'