Update PGS to contain tree based lithotype rules

examples/11_plurigaussian/07_tree_pgs.py

@@ -0,0 +1,133 @@
+"""
+PGS through decision trees
+--------------------------
+
+In typical PGS workflows, the lithotype is defined through a spatial rule. A
+more flexible approach can be taken such that the lithotype is represented as a
+decision tree. More specifically, this is given as a binary tree, where each
+node is a decision based on the values of the spatial random
+fields [Sektnan et al., 2024](https://doi.org/10.1007/s11004-024-10162-5). The
+leaf nodes are then assigned a discrete value which is given to the cell that
+is being evaluated. Here, a simple example is provided showing how to use the
+tree based approach in conducting plurigaussian simulation.
+
+As in the previous examples, we define the simulation domain and generate the
+necessary spatial random fields.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 2
+
+# no. of cells in both dimensions
+N = [150, 150]
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+
+model = gs.Gaussian(dim=dim, var=1, len_scale=10)
+srf = gs.SRF(model)
+field1 = srf.structured([x, y], seed=215253419)
+field2 = srf.structured([x, y], seed=192534221)
+
+###############################################################################
+# Decisions within the tree are carried out through user defined functions.
+# In this way, the lithotype is defined in a continuous space, not relying on
+# discretization. In this example we will use an ellipse as a decision boundary.
+# The function accepts a data dictionary, which contains the values of the
+# spatial random fields, and the parameters of the ellipse.
+
+
+def ellipse(data, key1, key2, c1, c2, s1, s2, angle=0):
+    x, y = data[key1] - c1, data[key2] - c2
+
+    if angle:
+        theta = np.deg2rad(angle)
+        c, s = np.cos(theta), np.sin(theta)
+        x, y = x * c + y * s, -x * s + y * c
+
+    return (x / s1) ** 2 + (y / s2) ** 2 <= 1
+
+
+###############################################################################
+# The decision tree is defined as a dictionary, where each node is a dictionary
+# itself. The root node is the first decision, which branches into two nodes,
+# one for each possible outcome. The leaf nodes are the final decisions, which
+# assign a discrete value to the given cell. The `func` key in each decision
+# node contains the function to be called, and the `args` key contains the
+# arguments to be passed to the function. These arguments must match the
+# parameters of the function. The `yes_branch` and `no_branch` keys contain the
+# names of the nodes to follow based on the outcome of the decision. `root`
+# must be given, but all other node names are arbitrary.
+#
+# Here, we define a simple decision tree with two phases. The first node
+# checks if the point is inside an ellipse, and if it is, it assigns the value
+# 1 to the lithotype. If it is not, it assigns the value 0. The parameters
+# of the ellipse are defined in the `args` dictionary. The keys `key1` and
+# `key2` refer to the spatial random fields, which are used to define the
+# ellipse. In the algorithm, the fields are indexed as `Z1`, `Z2`, etc.,
+# depending on the order in which they are passed to the PGS object. In this
+# case, `Z1` refers to `field1` and `Z2` refers to `field2`. The parameters
+# `c1`, `c2`, `s1`, `s2`, and `angle` define the center, scale, and rotation of
+# the ellipse, respectively.
+#
+# Lithotype decision functions operate in the latent Gaussian domain [-3,3]×[-3,3].
+# Since each GRF is N(0,1), approximately 99.7% of its values lie within ±3σ,
+# making ±3 a natural “full” range for defining splits or shapes.
+
+config = {
+    "root": {
+        "type": "decision",
+        "func": ellipse,
+        "args": {
+            "key1": "Z1",
+            "key2": "Z2",
+            "c1": 0,
+            "c2": 0,
+            "s1": 2.5,
+            "s2": 0.8,
+            "angle": -45,
+        },
+        "yes_branch": "phase1",
+        "no_branch": "phase0",
+    },
+    "phase1": {"type": "leaf", "action": 1},
+    "phase0": {"type": "leaf", "action": 0},
+}
+
+###############################################################################
+# With the tree configuration ready, we can create the PGS object as normal,
+# passing our domain size and spatial random fields. When generating the
+# plurigaussian field, we pass the tree configuration so GSTools knows which
+# PGS process to follow.
+
+pgs = gs.PGS(dim, [field1, field2])
+
+P = pgs(tree=config)
+
+###############################################################################
+# We can also compute the equivalent spatial lithotype.
+#
+# NB: If we want to compute L before P, we must pass the tree config to
+# `pgs.compute_lithotype(tree=config)`.
+
+L = pgs.compute_lithotype()
+
+###############################################################################
+# Finally, we plot `P` as well as the equivalent spatial lithotype `L`.
+
+fig, axs = plt.subplots(1, 2, figsize=(10, 5), dpi=150)
+
+im0 = axs[0].imshow(L, cmap="copper", origin="lower")
+axs[0].set_title("L")
+im1 = axs[1].imshow(P, cmap="copper", origin="lower")
+axs[1].set_title("P")
+
+###############################################################################
+#
+# .. image:: ../../pics/2d_tree_pgs.png
+#    :width: 400px
+#    :align: center

examples/11_plurigaussian/08_multi_field_pgs.py

@@ -0,0 +1,117 @@
+"""
+From bigaussian to plurigaussian simulation
+--------------------------------------------------
+
+In many PGS implementations, the dimensions of the simulation domain often
+matches the number of fields that are supplied. However, this is not a
+requirement of the PGS algorithm. In fact, it is possible to use multiple
+spatial random fields in PGS, which can be useful for more complex lithotype
+definitions. In this example, we will demonstrate how to use multiple SRFs in
+PGS. In GSTools, this is done by utlising the tree based architecture.
+
+Typically, PGS in two dimensions is carried out as a bigaussian simulation,
+where two random fields are used. Here, we will employ four. We begin by
+setting up the simulation domain and generating the necessary random fields,
+where the length scales of two of the fields are much larger than the other two.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 2
+
+N = [200, 200]
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+
+model = gs.Gaussian(dim=dim, var=1, len_scale=15)
+srf = gs.SRF(model)
+field1 = srf.structured([x, y], seed=201519)
+field2 = srf.structured([x, y], seed=199221)
+model = gs.Gaussian(dim=dim, var=1, len_scale=3)
+srf = gs.SRF(model)
+field3 = srf.structured([x, y], seed=1345134)
+field4 = srf.structured([x, y], seed=1351455)
+
+###############################################################################
+# As in the previous example, an ellipse is used as the decision boundary.
+
+
+def ellipse(data, key1, key2, c1, c2, s1, s2, angle=0):
+    x, y = data[key1] - c1, data[key2] - c2
+
+    if angle:
+        theta = np.deg2rad(angle)
+        c, s = np.cos(theta), np.sin(theta)
+        x, y = x * c + y * s, -x * s + y * c
+
+    return (x / s1) ** 2 + (y / s2) ** 2 <= 1
+
+
+###############################################################################
+# The configuration dictionary for the decision tree is defined as before, but
+# this time we pass the additional keys `Z3` and `Z4`, which refer to the
+# additional fields `field3` and `field4`. The decision tree is structured in a
+# way that the first decision node is based on the first two fields, and the
+# second decision node is based on the last two fields.
+
+config = {
+    "root": {
+        "type": "decision",
+        "func": ellipse,
+        "args": {
+            "key1": "Z1",
+            "key2": "Z2",
+            "c1": 0.7,
+            "c2": 0.7,
+            "s1": 1.3,
+            "s2": 1.3,
+        },
+        "yes_branch": "phase1",
+        "no_branch": "node1",
+    },
+    "node1": {
+        "type": "decision",
+        "func": ellipse,
+        "args": {
+            "key1": "Z3",
+            "key2": "Z4",
+            "c1": -0.7,
+            "c2": -0.7,
+            "s1": 1.3,
+            "s2": 1.3,
+            "angle": 30,
+        },
+        "yes_branch": "phase2",
+        "no_branch": "phase0",
+    },
+    "phase2": {"type": "leaf", "action": 2},
+    "phase1": {"type": "leaf", "action": 1},
+    "phase0": {"type": "leaf", "action": 0},
+}
+
+###############################################################################
+# Next, we create the PGS object, passing in all four fields, and generate the
+# plurigaussian field `P`
+
+pgs = gs.PGS(dim, [field1, field2, field3, field4])
+
+P = pgs(tree=config)
+
+###############################################################################
+# Finally, we plot `P`
+#
+# NB: In the current implementation, the calculation of the equivalent spatial
+# lithotype `L` is not supported for multiple fields
+
+plt.figure(figsize=(8, 6))
+plt.imshow(P, cmap="copper", origin="lower")
+
+###############################################################################
+#
+# .. image:: ../../pics/2d_multi_tree_pgs.png
+#    :width: 400px
+#    :align: center

examples/11_plurigaussian/09_3d_tree_pgs.py

@@ -0,0 +1,83 @@
+"""
+Three dimensional PGS through decision trees
+--------------------------------------------
+
+Let's apply the decision tree approach to three dimensional PGS
+"""
+
+import numpy as np
+
+import gstools as gs
+
+dim = 3
+
+# no. of cells in all dimensions
+N = [60] * dim
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+z = np.arange(N[2])
+
+###############################################################################
+# As we want to generate a three dimensional field, we generate three spatial
+# random fields (SRF) as our input. In this example, the number of fields
+# is equal to the number of dimensions, but as before, this is not a
+# requirement.
+
+model = gs.Gaussian(dim=dim, var=1, len_scale=[15, 15, 15])
+srf = gs.SRF(model)
+field1 = srf.structured([x, y, z], seed=20277519)
+field2 = srf.structured([x, y, z], seed=19727221)
+field3 = srf.structured([x, y, z], seed=21145612)
+
+###############################################################################
+# The decision tree will now utilise an ellipsoid as the decision boundary,
+# having a similar structure as the ellipse in the previous example.
+
+
+def ellipsoid(data, key1, key2, key3, c1, c2, c3, s1, s2, s3):
+    return ((data[key1] - c1) / s1) ** 2 + ((data[key2] - c2) / s2) ** 2 + (
+        (data[key3] - c3) / s3
+    ) ** 2 <= 1
+
+
+config = {
+    "root": {
+        "type": "decision",
+        "func": ellipsoid,
+        "args": {
+            "key1": "Z1",
+            "key2": "Z2",
+            "key3": "Z3",
+            "c1": 0,
+            "c2": 0,
+            "c3": 0,
+            "s1": 3,
+            "s2": 1,
+            "s3": 0.4,
+        },
+        "yes_branch": "phase1",
+        "no_branch": "phase0",
+    },
+    "phase0": {"type": "leaf", "action": 0},
+    "phase1": {"type": "leaf", "action": 1},
+}
+
+###############################################################################
+# As before, we initialise the PGS process, generate `P`, and plot to
+# visualise the results
+
+import pyvista as pv
+
+pgs = gs.PGS(dim, [field1, field2, field3])
+P = pgs(tree=config)
+
+grid = pv.ImageData(dimensions=N)
+grid.point_data["PGS"] = P.reshape(-1)
+grid.threshold(0.5, scalars="PGS").plot()
+
+###############################################################################
+#
+# .. image:: ../../pics/3d_tree_pgs.png
+#    :width: 400px
+#    :align: center

examples/11_plurigaussian/README.rst

@@ -1,14 +1,27 @@
-Plurigaussian Field Generation
-==============================
+Plurigaussian Simulation
+========================
 
-Plurigaussian field simulations (PGS) are used to simulate correlated fields
+Plurigaussian simulation (PGS) is used to simulate correlated fields
 of categorical data, e.g. lithofacies, hydrofacies, soil types, or
 cementitious materials.
-PGS uses one spatial random field (SRF) per dimension, e.g. two SRFs, when
-working with two dimensional data. Furthermore, PGS needs a field, which
-describes the categorical data and its spatial relations.
-This might sound more complicated then it is, as we will see in the following
-examples.
+In general, we define a categorical rule which dictates the relative
+frequency and connectivity of the phases present in the final realisation.
+We employ spatial random fields (SRFs) to map to this rule.
+This mapping determines the phase of a given point in the simulation domain.
+The definition of this rule limits how we can interact with it.
+For example, the rule may be defined spatially (e.g. as an image or array)
+or via a decision tree. Both forms will be explored in the following
+examples, highlighting their differences.
+Many PGS approaches constrain the number of input SRFs to match the
+dimensions of the simulation domain. This constraint stems from the
+implementation, not the method itself.
+With a spatial lithotype, we perform bigaussian and trigaussian
+simulations for two- and three-dimensional realisations, respectively.
+In contrast, the tree-based approach allows an arbitrary number of SRFs,
+yielding a fully *pluri*gaussian simulation.
+This may sound more complicated than it is; we will clarify everything
+in the examples that follow.
+
 
 Examples
 --------