Check for observartion point inside tesseroid eats too much memory

[santisoler]

Description of the problem:

The tesseroid forward modelling function checks if no observation point falls inside any tesseroid through the _check_points_outside_tesseroids function:

harmonica/harmonica/_forward/_tesseroid_utils.py

Lines 405 to 461 in b8f037e

	def _check_points_outside_tesseroids(coordinates, tesseroids):
	"""
	Check if computation points are not inside the tesseroids
	Parameters
	----------
	coordinates : 2d-array
	Array containing the coordinates of the computation points in the
	following order: ``longitude``, ``latitude`` and ``radius``.
	Both ``longitude`` and ``latitude`` must be in degrees.
	The array must have the following shape: (3, ``n_points``), where
	``n_points`` is the total number of computation points.
	tesseroids : 2d-array
	Array containing the boundaries of the tesseroids in the following
	order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
	Longitudinal and latitudinal boundaries must be in degrees.
	The array must have the following shape: (``n_tesseroids``, 6), where
	``n_tesseroids`` is the total number of tesseroids.
	This array of tesseroids must have longitude continuity and valid
	boundaries.
	Run ``_check_tesseroids`` before.
	"""
	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():
	err_msg = (
	"Found computation point inside tesseroid. "
	+ "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, :]
	)
	raise ValueError(err_msg)

This functions makes use of Numpy arrays to run these checks, but this is super inefficient, since Numpy will try to allocate massive arrays of shapes (n_coords, n_tesseroids).

Solution
We should rewrite this function using Numba instead, making sure that we aren't allocating any additional array and the function just checks if observation points don't fall inside any tesseroids. We could also benefit from Numba parallelizations and run the first for loop in parallel (this should follow the user's choice set through the parallel argument in tesseroid_gravity).

As a quick design, we could do something like:

@numba.jit(nopython=True, parallel=True)
def check_points_outside_tesseroids(coordinates, tesseroids):
    """..."""
    longitude, latitude, radius = coordinates[:]
    for i in numba.prange(longitude.size):
        for j in range(tesseroids.shape[0]):
            _check_point_inside_tesseroid(
                longitude[i],
                latitude[i],
                radius[i],
                tesseroids[j, 0],
                tesseroids[j, 1],
                tesseroids[j, 2],
                tesseroids[j, 3],
                tesseroids[j, 4],
                tesseroids[j, 5],
            )
            
@numba.jit(nopython=True)
def _check_point_inside_tesseroid(
    longitude, latitude, radius, west, east, south, north, bottom, top
):
    """Raise error if observation point falls inside tesseroid"""
    ...
    if ... :
        raise ValueError(...)

Check out https://github.com/fatiando/harmonica/blob/main/harmonica/_forward/prism.py for inspiration on how to code this.

Temporary solution

For the moment, the check can be bypassed by passing disable_checks=True, but it also disables other checks.