Skip to content

Commit 84a0010

Browse files
committed
[ENH] Add new likelihood method and corresponding test
Add `apparent_thickness_likelihood_II` function for calculating apparent thickness with structured grids. Introduce a new test case to validate functionality and enhance GemPy posterior visualization. Update documentation logs to reflect development progress.
1 parent a158d1b commit 84a0010

File tree

6 files changed

+118
-8
lines changed

6 files changed

+118
-8
lines changed

docs/dev_logs/2025_May.md

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,10 @@
1616
- [ ] Add Simon's example
1717

1818
### Doing:
19-
- [ ] Compressing code
20-
- [ ] GemPy posterior visualization
19+
- [x] Compressing code
20+
- [x] GemPy posterior visualization
21+
- [ ] Add proper thickness likelihood
22+
- [ ] Speed
2123

2224
### Possible Design:
2325
- Geological Model
Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,3 @@
1-
from _apparent_thickness import apparent_thickness_likelihood
2-
from _gravity_inv import gravity_inversion_likelihood
1+
from ._apparent_thickness import apparent_thickness_likelihood, apparent_thickness_likelihood_II
2+
from ._gravity_inv import gravity_inversion_likelihood
3+

gempy_probability/modules/likelihoods/_apparent_thickness.py

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,3 +15,14 @@ def apparent_thickness_likelihood(model_solutions: gp.data.Solutions) -> torch.T
1515
value=thickness.detach() # * This is only for az to track progress
1616
)
1717
return thickness
18+
19+
def apparent_thickness_likelihood_II(model_solutions: gp.data.Solutions, distance: float, element_lith_id: float) -> torch.Tensor:
20+
# TODO: element_lith_id should be an structured element
21+
simulated_well = model_solutions.octrees_output[0].last_output_center.custom_grid_values
22+
23+
# Create a boolean mask for all values between element_lith_id and element_lith_id+1
24+
mask = (simulated_well >= element_lith_id) & (simulated_well < (element_lith_id + 1))
25+
26+
# Count these values and compute thickness as the product of the count and the distance
27+
thickness = mask.sum() * distance
28+
return thickness

tests/test_likelihoods/__init__.py

Whitespace-only changes.
Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
import numpy as np
2+
import os
3+
4+
import gempy as gp
5+
import gempy_viewer as gpv
6+
from gempy.core.data.grid_modules import CustomGrid
7+
from gempy_engine.core.backend_tensor import BackendTensor
8+
from gempy_probability.modules.likelihoods import apparent_thickness_likelihood_II
9+
10+
11+
def test_gempy_posterior_I():
12+
current_dir = os.path.dirname(os.path.abspath(__file__))
13+
data_path = os.path.abspath(os.path.join(current_dir, '..', '..', 'examples', 'tutorials', 'data'))
14+
geo_model = gp.create_geomodel(
15+
project_name='Wells',
16+
extent=[0, 12000, -500, 500, 0, 4000],
17+
refinement=4,
18+
importer_helper=gp.data.ImporterHelper(
19+
path_to_orientations=os.path.join(data_path, "2-layers", "2-layers_orientations.csv"),
20+
path_to_surface_points=os.path.join(data_path, "2-layers", "2-layers_surface_points.csv")
21+
)
22+
)
23+
24+
BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH)
25+
26+
geo_model.interpolation_options.uni_degree = 0
27+
geo_model.interpolation_options.mesh_extraction = False
28+
29+
x_loc = 6000
30+
y_loc = 0
31+
z_loc = np.linspace(0, 4000, 100)
32+
dz = (4000 - 0)/100
33+
xyz_coord = np.array([[x_loc, y_loc, z] for z in z_loc])
34+
35+
custom_grid = CustomGrid(xyz_coord)
36+
geo_model.grid.custom_grid = custom_grid
37+
gp.set_active_grid(
38+
grid=geo_model.grid,
39+
grid_type=[gp.data.Grid.GridTypes.CUSTOM]
40+
)
41+
gp.compute_model(gempy_model=geo_model)
42+
p2d = gpv.plot_2d(
43+
model=geo_model,
44+
show_topography=False,
45+
legend=False,
46+
show_lith=False,
47+
show=False
48+
)
49+
50+
element_lith_id = 1.5
51+
thickness = apparent_thickness_likelihood_II(
52+
model_solutions=geo_model.solutions,
53+
distance=dz,
54+
element_lith_id=element_lith_id
55+
)
56+
57+
print(thickness)
58+
59+
p2d = gpv.plot_2d(geo_model, show_topography=False, legend=False, show=False)
60+
# --- New code: plot all custom grid points with lithology 1 ---
61+
#
62+
# Extract the simulated well values from the solution obtained.
63+
simulated_well = geo_model.solutions.octrees_output[0].last_output_center.custom_grid_values
64+
#
65+
# Create a boolean mask for all grid points whose value is between 1 and 2.
66+
mask = (simulated_well >= element_lith_id) & (simulated_well < (element_lith_id + 1))
67+
#
68+
# Retrieve the custom grid coordinates.
69+
grid_coords = geo_model.grid.custom_grid.values
70+
#
71+
# Use the mask to get only the grid points corresponding to lith id 1.
72+
# Note: Convert the boolean mask to a NumPy array if necessary.
73+
coords_lith1 = grid_coords[mask]
74+
#
75+
# Plot these points on the figure (using the x and z columns for a 2D image).
76+
p2d.axes[0].scatter(
77+
coords_lith1[:, 0],
78+
coords_lith1[:, 2],
79+
marker='o',
80+
s=80,
81+
c='green',
82+
label='Lith id = 1',
83+
zorder=9
84+
)
85+
p2d.axes[0].legend()
86+
# --------------------------------------------------------------
87+
88+
# Plot device location for reference.
89+
device_loc = np.array([[6e3, 0, 3700]])
90+
p2d.axes[0].scatter(
91+
device_loc[:, 0],
92+
device_loc[:, 2],
93+
marker='x',
94+
s=400,
95+
c='#DA8886',
96+
zorder=10
97+
)
98+
p2d.fig.show()

tests/test_posterior_analysis/test_gempy_posterior_I.py

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ def test_gempy_posterior_I():
2222
path_to_surface_points=os.path.join(data_path, "2-layers", "2-layers_surface_points.csv")
2323
)
2424
)
25-
25+
2626
BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH)
2727

2828
geo_model.interpolation_options.uni_degree = 0
@@ -47,7 +47,7 @@ def test_gempy_posterior_I():
4747
xyz[:, 2] = posterior_top_mean_z
4848
world_coord = geo_model.input_transform.apply_inverse(xyz)
4949
i = 0
50-
for i in range(0, 200, 5):
50+
for i in np.linspace(0, 200, 20).astype(int):
5151
gp.modify_surface_points(
5252
geo_model=geo_model,
5353
slice=0,
@@ -72,8 +72,6 @@ def test_gempy_posterior_I():
7272

7373
p2d.fig.show()
7474

75-
pass
76-
7775

7876
def _plot_posterior(data):
7977
az.plot_trace(data)

0 commit comments

Comments
 (0)