Add option to discard thin prisms of a layer when forward modelling

[mdtanker]

Add a new thickness_threshold parameter to the prism_layer.gravity() method that allows the user to ignore prisms below a certain thickness. This is implemented within a new _discard_thin_prisms(), function which is similar to the _discard_null_prisms() function. Add tests that check whether the calculated gravity is the same from manually removing prisms and removing prisms based on a threshold, and whether providing no threshold still produces the correct gravity.

I've added a parameter thickness_threshold to the prism_layer.gravity() function which allows the user to ignore prisms below a certain thickness. This is implemented within a new function _discard_thin_prisms(), which is similar to the function _discard_null_prisms() here.

I've also added a test that checks:

• whether the calculated gravity is the same from manually removing prisms and removing prisms based on a threshold

• whether providing no threshold still produces the correct gravity.

This shouldn't affect performance if no threshold is supplied since the discard is skipped if thickness_threshold is None.

I'm not sure what the performance hit is when a threshold is supplied.

Relevant issues/PRs:

Fixes #371

[mdtanker]

@santisoler thanks for reviewing so fast! Good catch with adding that additional test. I'm having trouble using the _discard_thin_prisms() function outside of prism_layer.py (i.e. inside the test_prism_layer.py). I've tried from .. import _discard_thin_prisms, which seems to work for prism_layer but not this function.

Also, I noticed most of the other functions in prism_layer are actually methods of the class DatasetAccessorPrismLayer. Would it make since to re-write this function as one of these methods?

[santisoler]

Hi @mdtanker!

I'm having trouble using the _discard_thin_prisms() function outside of prism_layer.py (i.e. inside the test_prism_layer.py). I've tried from .. import _discard_thin_prisms, which seems to work for prism_layer but not this function.

When you are trying to import with from .. import _discard_thin_prisms you are asking to go up one directory in the package (two dots .. mean "the parent directory" in Unix-like OSs and also in the import statements in Python) and import _discard_thin_prisms from there, but that function is not available at that level. Since our tree is:

├── harmonica
│   │ ...
│   ├── tests/
│   │   ├── test_prism_layer.py
│   │   ├── ...
│   ├── __init__.py
│   │ ...

when you run from .. import ... you are going to be able to import only the functions and classes that are listed in the harmonica/__init__.py. But since _discard_thin_prisms is a private function it won't be available from there.

The good thing is that you can always import private stuff, but you need to provide the whole route to them. In this case you can import the _discard_thin_prisms function with:

from ..forward.prism_layer import _discard_thin_prisms

Also, I noticed most of the other functions in prism_layer are actually methods of the class DatasetAccessorPrismLayer. Would it make since to re-write this function as one of these methods?

I'm glad you bring this question. While doing object-oriented designs is usually tempting to add a new method to the class, even if it's a private one. But this instinct can lead to bugs or even to repeat code that could be broadly used throughout the project in the future. Every time you need to decide if some function should be a method or a function ask yourself: is this function actually needing some attribute of the class? In this particular case the _discard_thin_prisms function only needs a set of prisms, their densities and the density threshold. It can be applied even to a general list of prisms: there's nothing that bounds this function to the prism layer class. So in my opinion this should be a private function for now.

By contrast, every other method of the prism layer is actually reading some attributes from self to return something new or modifying some subset of these attributes. That's why those are methods and not functions. But even those being methods is also an interface decision that could be discussed. It is better to run the forward of the prism layer with a method?

result = layer.prism_layer.gravity(coordinates, field)

or with a dedicated function?

result = prism_layer_gravity(coordinates, layer, field)

Probably not the discussion for this PR, but I just wanted to mention it since behind our designs we are always taking these types of decisions.

[mdtanker]

@santisoler thanks for the help. I've added a separate test to check that _discard_thin_prisms works correctly. Let me know if there's anything else you want to be tested here!

Thanks for the explanation on methods vs. functions! Very helpful to have that explained.

[mdtanker]

Thanks for those suggestions. I've added them, and a few other minor changes. Happy for you to merge whenever 👍

[santisoler]

Thanks @mdtanker! I think this is ready to be merged! 🚀