Forward modelling for point masses

[santisoler]

Add a forward modeling for point masses in geographic coordinates.
The coordinates of the computation points can be passed as list of arrays ([longitude, latitude, height]) while the point mass location can be specified as a list of its coordinates.
The user can specify which coordinate system is using for these coordinates (geodetic or spherical geocentric).
The time demanding portions of the code have been decorated with Numba's jit.
The fastmath argument is not used because it was proved to sacrifice accuracy over performance making test functions to fail.

Fixes #47

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.

[santisoler]

@leouieda I drafted portions of the point mass forward modelling code.
Regarding code style and user interface, what do you think?
I still need to make the code suitable for multiple point masses, but I want to have your opinion before getting too deep into the rabbit hole.

[leouieda]

Hi Santi, thanks for putting this together. I don't think we should be doing coordinate conversions in these functions. It's not much work for the user to convert coordinates before passing them in and makes sure they know what is going on. Plus, it's probably a bit wasteful to convert coordinates every time we call this function. Better to do it once outside of the function. This will probably be even more so when we do inversions.

Are you thinking of having this function calculate for multiple point masses? I'm not entirely sure how to handle this. We can probably optimize some things if we do this, like calculating sin and cos of the observation point latitude. But we also loose the simplicity of these functions. If the point masses have this, we should also do things this way for prisms, tesseroids, etc.

One advantage of leaving this as single element functions is that we can build other functions that compute on multiple elements and call these under the hood. For example, when calculating a regular grid of masses on a regular grid, there optimizations that we can make regarding symmetric points.

The main problem is that it's hard to do something generic that will work every case.

[santisoler]

Hi Santi, thanks for putting this together. I don't think we should be doing coordinate conversions in these functions. It's not much work for the user to convert coordinates before passing them in and makes sure they know what is going on. Plus, it's probably a bit wasteful to convert coordinates every time we call this function. Better to do it once outside of the function. This will probably be even more so when we do inversions.

You're right, as I pointed on #46. I'll remove the coordinate conversion and assume spherical coordinates.

[santisoler]

Are you thinking of having this function calculate for multiple point masses? I'm not entirely sure how to handle this. We can probably optimize some things if we do this, like calculating sin and cos of the observation point latitude. But we also loose the simplicity of these functions. If the point masses have this, we should also do things this way for prisms, tesseroids, etc.

One advantage of leaving this as single element functions is that we can build other functions that compute on multiple elements and call these under the hood. For example, when calculating a regular grid of masses on a regular grid, there optimizations that we can make regarding symmetric points.

The main problem is that it's hard to do something generic that will work every case.

When we build forward models we have two summations: one over the sources and one over the computation points. For example:

for source in sources:
    for point in computation_points:
        result += forward_model(computation_point, source)

If we leave this functions with single element computation we are speeding up the nested for through Numba, but the first one will be a slow Python loop. That's why I always thought we should implement both summations on fast execution code.

Did you consider this? Or you were thinking to jit the point_mass_gravity() function as well?

Regarding the multiple elements implementation, I don't think it would be so hard. Instead of asking for a single point mass, the function could take a collection of them (in the same way the coordinates are passed), ravel the point masses coordinates and pass them to the jitted function. In there we just need to add an outer loop that iterates over each point mass. Although we are adding a few lines, I don't think this would make the function loose simplicity.

P.S: In my opinion, the forward models that we'll build should have the same logic, i.e. if one works with multiple elements, they all should also.

[santisoler]

@leouieda Can we continue discussing this PR so I can keep writing with a more precise North?

[santisoler]

When trying to get the test coverage, pytest cannot get into the jitted functions.
It would be nice to run tests on both jitted and not jitted versions of the functions.
Disabling jit will allow pytest to get proper coverage, while enabling it will guarantee that the functions that ultimately users will use work as expected (for example, setting fastmath=True sacrifices accuracy, and it's important to test if that sacrifice is significant).

According to Numba documentation this can be done setting the environment variable NUMBA_DISABLE_JIT to 1, although I haven't been able to change its value inside the test functions in order to run them twice: one with NUMBA_DISABLE_JIT=0 and NUMBA_DISABLE_JIT=1.
A simple solution would be to setting this environment variable on the Makefile and run pytest twice, but I don't like it very much.

[santisoler]

The implementation of this enhancement has been moved to #60