Rewrite check for point inside tesseroid with Numba (#419)

Rewrite the functions that checks if a computation point falls inside a
tesseroid using Numba instad of Numpy. This greatly reduce the memory
need of the Numpy functions, that required to build large matrices even
for medium-size problems.

harmonica/_forward/_tesseroid_utils.py

@@ -402,7 +402,7 @@ def _check_tesseroids(tesseroids):
     return tesseroids
 
 
-def _check_points_outside_tesseroids(coordinates, tesseroids):
+def check_points_outside_tesseroids(coordinates, tesseroids):
     """
     Check if computation points are not inside the tesseroids
 
@@ -423,44 +423,53 @@ def _check_points_outside_tesseroids(coordinates, tesseroids):
         This array of tesseroids must have longitude continuity and valid
         boundaries.
         Run ``_check_tesseroids`` before.
+
+    Raises
+    ------
+    ValueError
+        If any computation point falls inside any tesseroid.
     """
-    longitude, latitude, radius = coordinates[:]
-    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
-    # Longitudinal boundaries of the tesseroid must be compared with
-    # longitudinal coordinates of computation points when moved to
-    # [0, 360) and [-180, 180).
-    longitude_360 = longitude % 360
-    longitude_180 = ((longitude + 180) % 360) - 180
-    inside_longitude = np.logical_or(
-        np.logical_and(
-            west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east
-        ),
-        np.logical_and(
-            west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east
-        ),
-    )
-    inside_latitude = np.logical_and(
-        south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north
-    )
-    inside_radius = np.logical_and(
-        bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top
-    )
-    # Build array of booleans.
-    # The (i, j) element is True if the computation point i is inside the
-    # tesseroid j.
-    inside = inside_longitude * inside_latitude * inside_radius
-    if inside.any():
+    longitude, latitude, radius = coordinates
+    conflicting = _check_points_outside_tesseroids(coordinates, tesseroids)
+    if conflicting:
         err_msg = (
-            "Found computation point inside tesseroid. "
-            + "Computation points must be outside of tesseroids.\n"
+            "Found computation point(s) inside tesseroid(s). "
+            "Computation points must be outside of tesseroids.\n"
         )
-        for point_i, tess_i in np.argwhere(inside):
-            err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format(
-                coordinates[:, point_i], tesseroids[tess_i, :]
+        for i, j in conflicting:
+            west, east, south, north, bottom, top = tesseroids[j, :]
+            err_msg += (
+                f" - Computation point '({longitude[i]}, {latitude[i]}, {radius[i]})' "
+                "inside tesseroid "
+                f"'({west}, {east}, {south}, {north}, {bottom}, {top})'.\n"
             )
         raise ValueError(err_msg)
 
 
+@jit(nopython=True)
+def _check_points_outside_tesseroids(coordinates, tesseroids):
+    """
+    Check if observation points fall inside tesseroids.
+    """
+    longitude, latitude, radius = coordinates[:]
+    conflicting = []
+    for i in range(longitude.size):
+        for j in range(tesseroids.shape[0]):
+            # Longitudinal boundaries of the tesseroid must be compared with
+            # longitudinal coordinates of computation points when moved to
+            # [0, 360) and [-180, 180).
+            longitude_360 = longitude[i] % 360
+            longitude_180 = ((longitude[i] + 180) % 360) - 180
+            west, east, south, north, bottom, top = tesseroids[j, :]
+            if (
+                (west < longitude_180 < east or west < longitude_360 < east)
+                and south < latitude[i] < north
+                and bottom < radius[i] < top
+            ):
+                conflicting.append((i, j))
+    return conflicting
+
+
 def _longitude_continuity(tesseroids):
     """
     Modify longitudinal boundaries of tesseroids to ensure longitude continuity

harmonica/_forward/tesseroid.py

@@ -13,9 +13,9 @@
 from ..constants import GRAVITATIONAL_CONST
 from ._tesseroid_utils import (
     _adaptive_discretization,
-    _check_points_outside_tesseroids,
     _check_tesseroids,
     _discard_null_tesseroids,
+    check_points_outside_tesseroids,
     gauss_legendre_quadrature,
     glq_nodes_weights,
 )
@@ -163,7 +163,7 @@ def tesseroid_gravity(
     # Sanity checks for tesseroids and computation points
     if not disable_checks:
         tesseroids = _check_tesseroids(tesseroids)
-        _check_points_outside_tesseroids(coordinates, tesseroids)
+        check_points_outside_tesseroids(coordinates, tesseroids)
     # Check if density is a function or constant values
     if callable(density):
         # Run density-based discretization on each tesseroid

harmonica/tests/test_tesseroid.py

@@ -24,8 +24,8 @@
     MAX_DISCRETIZATIONS,
     STACK_SIZE,
     _adaptive_discretization,
-    _check_points_outside_tesseroids,
     _check_tesseroids,
+    check_points_outside_tesseroids,
     tesseroid_gravity,
 )
 from ..constants import GRAVITATIONAL_CONST
@@ -129,7 +129,7 @@ def test_valid_tesseroid():
     _check_tesseroids(np.atleast_2d([w, e, s, n, bottom, bottom]))
     # Check if valid tesseroid with west > east is not caught
     _check_tesseroids(np.atleast_2d([350, 10, s, n, bottom, top]))
-    # Check if tesseroids on [-180, 180) are not caught
+    # Check if tesseroids on -180, 180 are not caught
     _check_tesseroids(np.atleast_2d([-70, -60, s, n, bottom, top]))
     _check_tesseroids(np.atleast_2d([-10, 10, s, n, bottom, top]))
     _check_tesseroids(np.atleast_2d([-150, 150, s, n, bottom, top]))
@@ -207,35 +207,71 @@ def test_disable_checks():
     npt.assert_allclose(invalid_result, -valid_result)
 
 
-def test_point_inside_tesseroid():
-    "Check if a computation point inside the tesseroid is caught"
-    tesseroids = np.atleast_2d([-10, 10, -10, 10, 100, 200])
-    # Test if outside point is not caught
-    points = [
-        np.atleast_2d([0, 0, 250]).T,  # outside point on radius
-        np.atleast_2d([20, 0, 150]).T,  # outside point on longitude
-        np.atleast_2d([0, 20, 150]).T,  # outside point on latitude
-        np.atleast_2d([0, 0, 200]).T,  # point on top surface
-        np.atleast_2d([0, 0, 100]).T,  # point on bottom surface
-        np.atleast_2d([-10, 0, 150]).T,  # point on western surface
-        np.atleast_2d([10, 0, 150]).T,  # point on eastern surface
-        np.atleast_2d([0, -10, 150]).T,  # point on southern surface
-        np.atleast_2d([0, 10, 150]).T,  # point on northern surface
-    ]
-    for coordinates in points:
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    # Test if computation point is inside the tesseroid
-    coordinates = np.atleast_2d([0, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    # Test if computation point with phased longitude is inside the tesseroid
-    coordinates = np.atleast_2d([360, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    tesseroids = np.atleast_2d([260, 280, -10, 10, 100, 200])
-    coordinates = np.atleast_2d([-90, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
+class TestPointInsideTesseroid:
+    """Test checks for points inside tesseroids."""
+
+    @pytest.mark.use_numba
+    def test_points_outside_tesseroid(self):
+        """Check if error is not raised when no point is inside a tesseroid."""
+        points = np.array(
+            [
+                [0, 0, 250],  # outside point on radius
+                [20, 0, 150],  # outside point on longitude
+                [0, 20, 150],  # outside point on latitude
+                [0, 0, 200],  # point on top surface
+                [0, 0, 100],  # point on bottom surface
+                [-10, 0, 150],  # point on western surface
+                [10, 0, 150],  # point on eastern surface
+                [0, -10, 150],  # point on southern surface
+                [0, 10, 150],  # point on northern surface
+            ]
+        ).T
+        tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200])
+        check_points_outside_tesseroids(points, tesseroid)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("points", [[0, 0, 150], [360, 0, 150]])
+    def test_point_inside_tesseroid(self, points):
+        """Check if error is raised when point fall inside the tesseroid."""
+        points = np.atleast_2d(points).T
+        tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200])
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(points, tesseroid)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize(
+        "point, tesseroid",
+        [
+            ([360, 0, 150], [-10, 10, -10, 10, 100, 200]),
+            ([-90, 0, 150], [260, 280, -10, 10, 100, 200]),
+        ],
+    )
+    def test_point_phased_longitude(self, point, tesseroid):
+        """Test if error is raised when the longitude coords are phased."""
+        point = np.atleast_2d(point).T
+        tesseroid = np.atleast_2d(tesseroid)
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(point, tesseroid)
+
+    @pytest.mark.use_numba
+    def test_multiple_points_and_tesseroids(self):
+        """Check if error is raised with multiple points and tesseroids."""
+        tesseroids = np.atleast_2d(
+            [
+                [-10, 10, -10, 10, 100, 200],
+                [20, 30, 20, 30, 400, 500],
+                [-50, -40, -30, -20, 100, 500],
+            ]
+        )
+        points = np.array(
+            [
+                [0, 0, 150],
+                [80, 82, 4000],
+                [10, 10, 450],
+            ]
+        ).T
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(points, tesseroids)
 
 
 @pytest.mark.use_numba