From fd8dc51ad2267ab34029bf2b5fb5ea7d902a1c99 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 4 May 2025 08:31:43 +0100 Subject: [PATCH 1/3] [WIP] Moving likelihood function to module --- gempy_probability/modules/likelihoods/__init__.py | 1 + .../modules/likelihoods/_apparent_thickness.py | 15 +++++++++++++++ .../modules/model_definition/model_examples.py | 14 ++------------ 3 files changed, 18 insertions(+), 12 deletions(-) create mode 100644 gempy_probability/modules/likelihoods/__init__.py create mode 100644 gempy_probability/modules/likelihoods/_apparent_thickness.py diff --git a/gempy_probability/modules/likelihoods/__init__.py b/gempy_probability/modules/likelihoods/__init__.py new file mode 100644 index 0000000..946c8d7 --- /dev/null +++ b/gempy_probability/modules/likelihoods/__init__.py @@ -0,0 +1 @@ +from _apparent_thickness import apparent_thickness_likelihood \ No newline at end of file diff --git a/gempy_probability/modules/likelihoods/_apparent_thickness.py b/gempy_probability/modules/likelihoods/_apparent_thickness.py new file mode 100644 index 0000000..12f4f16 --- /dev/null +++ b/gempy_probability/modules/likelihoods/_apparent_thickness.py @@ -0,0 +1,15 @@ +import gempy as gp +import torch +import pyro + +def apparent_thickness_likelihood(model_solutions: gp.data.Solutions) -> torch.Tensor: + """ + This function computes the thickness of the geological layer. + """ + simulated_well = model_solutions.octrees_output[0].last_output_center.custom_grid_values + thickness = simulated_well.sum() + pyro.deterministic( + name=r'$\mu_{thickness}$', + value=thickness.detach() # * This is only for az to track progress + ) + return thickness diff --git a/gempy_probability/modules/model_definition/model_examples.py b/gempy_probability/modules/model_definition/model_examples.py index 5765d86..fc8aeea 100644 --- a/gempy_probability/modules/model_definition/model_examples.py +++ b/gempy_probability/modules/model_definition/model_examples.py @@ -2,6 +2,7 @@ import gempy.core.data import gempy_engine from gempy.modules.data_manipulation.engine_factory import interpolation_input_from_structural_frame +from gempy_probability.modules.likelihoods import apparent_thickness_likelihood import pyro import pyro.distributions as dist @@ -53,7 +54,7 @@ def model(geo_model: gempy.core.data.GeoModel, normal, y_obs_list): # Compute and observe the thickness of the geological layer model_solutions: gp.data.Solutions = geo_model.solutions - thickness = _likelihood(model_solutions) + thickness = apparent_thickness_likelihood(model_solutions) # endregion @@ -67,14 +68,3 @@ def model(geo_model: gempy.core.data.GeoModel, normal, y_obs_list): ) -def _likelihood(model_solutions: gp.data.Solutions) -> torch.Tensor: - """ - This function computes the thickness of the geological layer. - """ - simulated_well = model_solutions.octrees_output[0].last_output_center.custom_grid_values - thickness = simulated_well.sum() - pyro.deterministic( - name=r'$\mu_{thickness}$', - value=thickness.detach() # * This is only for az to track progress - ) - return thickness From a158d1b9ae4aa042194994db62b2f453db1b56ea Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 4 May 2025 08:35:27 +0100 Subject: [PATCH 2/3] [WIP] Adding basic implementation of gravity --- gempy_probability/modules/likelihoods/__init__.py | 3 ++- .../modules/likelihoods/_apparent_thickness.py | 2 ++ gempy_probability/modules/likelihoods/_gravity_inv.py | 9 +++++++++ 3 files changed, 13 insertions(+), 1 deletion(-) create mode 100644 gempy_probability/modules/likelihoods/_gravity_inv.py diff --git a/gempy_probability/modules/likelihoods/__init__.py b/gempy_probability/modules/likelihoods/__init__.py index 946c8d7..a496a5e 100644 --- a/gempy_probability/modules/likelihoods/__init__.py +++ b/gempy_probability/modules/likelihoods/__init__.py @@ -1 +1,2 @@ -from _apparent_thickness import apparent_thickness_likelihood \ No newline at end of file +from _apparent_thickness import apparent_thickness_likelihood +from _gravity_inv import gravity_inversion_likelihood diff --git a/gempy_probability/modules/likelihoods/_apparent_thickness.py b/gempy_probability/modules/likelihoods/_apparent_thickness.py index 12f4f16..74cab7a 100644 --- a/gempy_probability/modules/likelihoods/_apparent_thickness.py +++ b/gempy_probability/modules/likelihoods/_apparent_thickness.py @@ -5,6 +5,8 @@ def apparent_thickness_likelihood(model_solutions: gp.data.Solutions) -> torch.Tensor: """ This function computes the thickness of the geological layer. + + Notes: This is not completed """ simulated_well = model_solutions.octrees_output[0].last_output_center.custom_grid_values thickness = simulated_well.sum() diff --git a/gempy_probability/modules/likelihoods/_gravity_inv.py b/gempy_probability/modules/likelihoods/_gravity_inv.py new file mode 100644 index 0000000..889af2d --- /dev/null +++ b/gempy_probability/modules/likelihoods/_gravity_inv.py @@ -0,0 +1,9 @@ +import pyro +import torch +import gempy as gp + +def gravity_inversion_likelihood(model_solutions: gp.data.Solutions) -> torch.Tensor: + simulated_geophysics = model_solutions.gravity + pyro.deterministic(r'$\mu_{gravity}$', simulated_geophysics) + + return simulated_geophysics \ No newline at end of file From 84a00100da0fa32b487767ec23dbe6060aa62955 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 4 May 2025 09:58:24 +0100 Subject: [PATCH 3/3] [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. --- docs/dev_logs/2025_May.md | 6 +- .../modules/likelihoods/__init__.py | 5 +- .../likelihoods/_apparent_thickness.py | 11 +++ tests/test_likelihoods/__init__.py | 0 .../test_geometric_likelihoods.py | 98 +++++++++++++++++++ .../test_gempy_posterior_I.py | 6 +- 6 files changed, 118 insertions(+), 8 deletions(-) create mode 100644 tests/test_likelihoods/__init__.py create mode 100644 tests/test_likelihoods/test_geometric_likelihoods.py diff --git a/docs/dev_logs/2025_May.md b/docs/dev_logs/2025_May.md index eb52903..57859b3 100644 --- a/docs/dev_logs/2025_May.md +++ b/docs/dev_logs/2025_May.md @@ -16,8 +16,10 @@ - [ ] Add Simon's example ### Doing: -- [ ] Compressing code -- [ ] GemPy posterior visualization +- [x] Compressing code +- [x] GemPy posterior visualization +- [ ] Add proper thickness likelihood +- [ ] Speed ### Possible Design: - Geological Model diff --git a/gempy_probability/modules/likelihoods/__init__.py b/gempy_probability/modules/likelihoods/__init__.py index a496a5e..e88d21a 100644 --- a/gempy_probability/modules/likelihoods/__init__.py +++ b/gempy_probability/modules/likelihoods/__init__.py @@ -1,2 +1,3 @@ -from _apparent_thickness import apparent_thickness_likelihood -from _gravity_inv import gravity_inversion_likelihood +from ._apparent_thickness import apparent_thickness_likelihood, apparent_thickness_likelihood_II +from ._gravity_inv import gravity_inversion_likelihood + diff --git a/gempy_probability/modules/likelihoods/_apparent_thickness.py b/gempy_probability/modules/likelihoods/_apparent_thickness.py index 74cab7a..40c6950 100644 --- a/gempy_probability/modules/likelihoods/_apparent_thickness.py +++ b/gempy_probability/modules/likelihoods/_apparent_thickness.py @@ -15,3 +15,14 @@ def apparent_thickness_likelihood(model_solutions: gp.data.Solutions) -> torch.T value=thickness.detach() # * This is only for az to track progress ) return thickness + +def apparent_thickness_likelihood_II(model_solutions: gp.data.Solutions, distance: float, element_lith_id: float) -> torch.Tensor: + # TODO: element_lith_id should be an structured element + simulated_well = model_solutions.octrees_output[0].last_output_center.custom_grid_values + + # Create a boolean mask for all values between element_lith_id and element_lith_id+1 + mask = (simulated_well >= element_lith_id) & (simulated_well < (element_lith_id + 1)) + + # Count these values and compute thickness as the product of the count and the distance + thickness = mask.sum() * distance + return thickness diff --git a/tests/test_likelihoods/__init__.py b/tests/test_likelihoods/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_likelihoods/test_geometric_likelihoods.py b/tests/test_likelihoods/test_geometric_likelihoods.py new file mode 100644 index 0000000..d6ee9ab --- /dev/null +++ b/tests/test_likelihoods/test_geometric_likelihoods.py @@ -0,0 +1,98 @@ +import numpy as np +import os + +import gempy as gp +import gempy_viewer as gpv +from gempy.core.data.grid_modules import CustomGrid +from gempy_engine.core.backend_tensor import BackendTensor +from gempy_probability.modules.likelihoods import apparent_thickness_likelihood_II + + +def test_gempy_posterior_I(): + current_dir = os.path.dirname(os.path.abspath(__file__)) + data_path = os.path.abspath(os.path.join(current_dir, '..', '..', 'examples', 'tutorials', 'data')) + geo_model = gp.create_geomodel( + project_name='Wells', + extent=[0, 12000, -500, 500, 0, 4000], + refinement=4, + importer_helper=gp.data.ImporterHelper( + path_to_orientations=os.path.join(data_path, "2-layers", "2-layers_orientations.csv"), + path_to_surface_points=os.path.join(data_path, "2-layers", "2-layers_surface_points.csv") + ) + ) + + BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH) + + geo_model.interpolation_options.uni_degree = 0 + geo_model.interpolation_options.mesh_extraction = False + + x_loc = 6000 + y_loc = 0 + z_loc = np.linspace(0, 4000, 100) + dz = (4000 - 0)/100 + xyz_coord = np.array([[x_loc, y_loc, z] for z in z_loc]) + + custom_grid = CustomGrid(xyz_coord) + geo_model.grid.custom_grid = custom_grid + gp.set_active_grid( + grid=geo_model.grid, + grid_type=[gp.data.Grid.GridTypes.CUSTOM] + ) + gp.compute_model(gempy_model=geo_model) + p2d = gpv.plot_2d( + model=geo_model, + show_topography=False, + legend=False, + show_lith=False, + show=False + ) + + element_lith_id = 1.5 + thickness = apparent_thickness_likelihood_II( + model_solutions=geo_model.solutions, + distance=dz, + element_lith_id=element_lith_id + ) + + print(thickness) + + p2d = gpv.plot_2d(geo_model, show_topography=False, legend=False, show=False) + # --- New code: plot all custom grid points with lithology 1 --- + # + # Extract the simulated well values from the solution obtained. + simulated_well = geo_model.solutions.octrees_output[0].last_output_center.custom_grid_values + # + # Create a boolean mask for all grid points whose value is between 1 and 2. + mask = (simulated_well >= element_lith_id) & (simulated_well < (element_lith_id + 1)) + # + # Retrieve the custom grid coordinates. + grid_coords = geo_model.grid.custom_grid.values + # + # Use the mask to get only the grid points corresponding to lith id 1. + # Note: Convert the boolean mask to a NumPy array if necessary. + coords_lith1 = grid_coords[mask] + # + # Plot these points on the figure (using the x and z columns for a 2D image). + p2d.axes[0].scatter( + coords_lith1[:, 0], + coords_lith1[:, 2], + marker='o', + s=80, + c='green', + label='Lith id = 1', + zorder=9 + ) + p2d.axes[0].legend() + # -------------------------------------------------------------- + + # Plot device location for reference. + device_loc = np.array([[6e3, 0, 3700]]) + p2d.axes[0].scatter( + device_loc[:, 0], + device_loc[:, 2], + marker='x', + s=400, + c='#DA8886', + zorder=10 + ) + p2d.fig.show() diff --git a/tests/test_posterior_analysis/test_gempy_posterior_I.py b/tests/test_posterior_analysis/test_gempy_posterior_I.py index 09e3a31..958b800 100644 --- a/tests/test_posterior_analysis/test_gempy_posterior_I.py +++ b/tests/test_posterior_analysis/test_gempy_posterior_I.py @@ -22,7 +22,7 @@ def test_gempy_posterior_I(): path_to_surface_points=os.path.join(data_path, "2-layers", "2-layers_surface_points.csv") ) ) - + BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH) geo_model.interpolation_options.uni_degree = 0 @@ -47,7 +47,7 @@ def test_gempy_posterior_I(): xyz[:, 2] = posterior_top_mean_z world_coord = geo_model.input_transform.apply_inverse(xyz) i = 0 - for i in range(0, 200, 5): + for i in np.linspace(0, 200, 20).astype(int): gp.modify_surface_points( geo_model=geo_model, slice=0, @@ -72,8 +72,6 @@ def test_gempy_posterior_I(): p2d.fig.show() - pass - def _plot_posterior(data): az.plot_trace(data)