Skip to content

Plurigaussian fields #370

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 46 commits into from
Apr 27, 2025
Merged

Plurigaussian fields #370

merged 46 commits into from
Apr 27, 2025

Conversation

LSchueler
Copy link
Member

Plurigaussian simulations (PGS) are a great way to easily increase the flexibility of Gaussian random fields. With this PR I want to directly incorporate them into GSTools. It not only includes the implementation, but also unittests and a few examples, with which users can familiarize themselves with PGS, as I do not find them very intuitive at first.

However, I still have a few open questions:

  • Should PGS really be a class? I mostly did this for GSTools to keep a more or less constant structure. But I don't see any benefits of using a class here. It only makes it a bit more complicated to use PGSs, as you have to first create an instance and then call it, instead of only calling a pgs function. Any opinion on this point @MuellerSeb ?

  • At one point we might have to think about a new structure for the examples. In 01_random_field we have examples for general field generation, but based on the randomization method and some examples showing the Fourier method. Then we have 04_vector_field with the vector field generation. And now we have 11_plurigaussian. What do you think @MuellerSeb , should we discuss this soon?

  • Is the argument facies the best name for a method which is used in many different fields?

  • Any other paper(s) we should cite?

@LSchueler LSchueler added enhancement New feature or request help wanted Extra attention is needed labels Nov 19, 2024
@LSchueler LSchueler added this to the 1.7 milestone Nov 19, 2024
@LSchueler LSchueler requested a review from MuellerSeb November 19, 2024 17:20
@LSchueler LSchueler self-assigned this Nov 19, 2024
I'm not able to install old numpy versions locally and the exception
documentation of numpy is non existent and I can't bother to go through
its source code, so I'll just see what the actions show.
I don't understand why `np.AxisError` doesn't work for np 2.x
@LSchueler
Copy link
Member Author

LSchueler commented Nov 20, 2024

Wow, that's a curious detail I stumbled upon...

And for some reason, the old numpy.AxisError didn't work for me with numpy 2.x, which was moved to numpy.exceptions.AxisError, but for backwards compatibility was supposed to be kept in the old place too.

@LSchueler
Copy link
Member Author

Thanks to @EJRicketts's feedback, I updated the examples and fixed a bug, where the L field could have an offset.

@EJRicketts
Copy link

Thank you or sharing this with me, it's nice to see that the offset has been sorted :) This was my only major comment.

In general, I think it would be nice to include an example that shows the use of conditional random fields and also periodic random fields.

For field scale problems with associated experimental observations, combining a well chosen lithotype with conditional random fields can be quite powerful in getting nice realistic representations.

With respect to the periodic case, periodicity in the input fields will result in periodicity in the final field, which for my research purposes have been useful in material characterisation. I published some work on this for cementitious materials recently: https://doi.org/10.1007/s11242-024-02074-z
It would be easy to replicate these results for the example, but perhaps a more geoenvironmental example would be best! In any case, a small mention of the paper in the example would be much appreciated 😊

Regarding citing works, there is a nice book on PGS: https://link.springer.com/book/10.1007/978-3-642-19607-2
The authors are (mostly) from the birth place of PGS (Mines-Paristech), some of which are the original authors who took the idea of truncated random fields and extended it to PGS. It feels important to include this.
Xavier Emery has also had quite a big influence on PGS, so it might be worth citing their work also: https://doi.org/10.1016/j.cageo.2007.01.006

@LSchueler
Copy link
Member Author

@MuellerSeb This is ready for review

Copy link
Member

@MuellerSeb MuellerSeb left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is awesome. And so cool to see what can be done in geo-statistics with some sugar on it.

A well done contribution and beautiful examples. 🎉
I would only like to see better naming for the used variables and especially for method arguments and attributes. Maybe this is a common thing in the pgs community, but I would like to have descriptive names in GSTools.

What do you think @LSchueler ?

Comment on lines 80 to 82
# .. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/3d_pgs.png
# :width: 400px
# :align: center
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not shown in the docs.

Copy link
Member Author

@LSchueler LSchueler Apr 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think, I hope this is shown as soon as this branch is merged into main. At least that's how I remember I did it the last time, I uploaded an image to the docs.
Can you confirm this?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, we shouldn't use html adresses for images stored in docs folder (as you can see by the use of "master"). We did that in the past due to stupidity. Sphinx should be capable of using local paths to reference images.

Comment on lines 39 to 83
# no. of grid cells of L-field
M = [60, 50]

# size of the rectangles
R = [10, 8]

# positions of some of the shapes for concise indexing
S1 = [1, -9]
S2 = [-5, 3]
S3 = [-5, -5]

L = np.zeros(M)
# a small upper triangular helper matrix to create the triangle
tri = np.triu(np.ones((R[0], R[0])))
# the triangle
L[
M[0] // 2 + S1[0] : M[0] // 2 + S1[0] + R[0],
M[1] // 2 + S1[1] : M[1] // 2 + S1[1] + R[0],
] = tri
# the first rectangle
L[
M[0] // 2 + S2[0] - R[0] // 2 : M[0] // 2 + S2[0] + R[0] // 2,
M[1] // 2 + S2[1] - R[1] // 2 : M[1] // 2 + S2[1] + R[1] // 2,
] = 2
# the second rectangle
L[
M[0] // 2 + S3[0] - R[0] // 2 : M[0] // 2 + S3[0] + R[0] // 2,
M[1] // 2 + S3[1] - R[1] // 2 : M[1] // 2 + S3[1] + R[1] // 2,
] = 3
# some very narrow rectangles
for i in range(4):
L[
M[0] // 2 + S1[0] : M[0] // 2 + S1[0] + R[0],
M[1] // 2
+ S1[1]
+ R[1]
+ 3
+ 2 * i : M[1] // 2
+ S1[1]
+ R[1]
+ 4
+ 2 * i,
] = (
4 + i
)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know if this is common in the literature about PGS, but I really dislike all the single uppercase Letters for matrices used in this context.

I think the whole API and documentation of GSTools is quiet descriptive, or at least tries to be, and this reads almost like Fortran 77 code.

Couldn't we come up with something like a descriptive and self explaining naming convention for these "things"?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ouch, that Fortran 77 comparison hurt my feelings :-D

I found that code pretty ugly already and didn't want to reduce readability by introducing longer variable names. But changed most of the variable names to longer and more descriptive ones, because your criticism is absolutely valid. I didn't change all of them and kept e.g. M, similarly as we often use N for the number of grid points.
Are you okay with this compromise?

@LSchueler
Copy link
Member Author

If you think that the image will be shown once this is merged and you like my new variable names, this could be merged?

Copy link
Member

@MuellerSeb MuellerSeb left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fixed the image issue in the documentation.

So now: LGTM

In a next step, we should provide a field child-class using PGS internally. For now this is a cool bare-bone feature.

Thanks for this neat work!

@LSchueler LSchueler merged commit 7da6000 into main Apr 27, 2025
32 of 33 checks passed
@LSchueler LSchueler deleted the pgs branch April 27, 2025 19:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request help wanted Extra attention is needed
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants