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

[santisoler]

Thanks @mdtanker for opening this PR! This is looking great!

I left a few little comments below, nothing too big. I think the only thing that we would like to improve before merging this are the tests. We are currently testing if the gravity fields of two layers (one that have discarded thin prisms and one that doesn't) generate the same gravity field. This test makes sense, but we are not actually testing if discarding prisms based on thickness actually works: assume that the function has a bug and it doesn't discard any prism, the gravity field would still be exactly the same and all tests will pass, but the bug will still be there.

Would you like to write an extra test that checks just the private _discard_thin_prisms function? Basically it should build a set of prisms with different thicknesses and a set of densities (would be better if they all have different densities), pass them to the _discard_thin_prisms function and check if the output doesn't contain the thin prisms and the density values are the correct ones: the ones that correspond with the thick prisms.

The advantage of writing private functions in a modular way is that it's way easier to test them independently of the higher level code that uses them. This is usually known as unit testing where we test the smallest units of our codebase independently of the rest.

[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.

[santisoler]

Thanks @mdtanker for adding that test! I left a few comments below, mainly about the docstrings. After those are fixed I think this should be ready to be merged.

[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! 🚀