[DOC] Getting the examples running

Author: Leguark

docs/dev_logs/2025_May.md

@@ -44,4 +44,11 @@
 
 ### Likelihood functions
 - ? Should I start to create a catalog of likelihood functions that goes with gempy?
-  - This would be useful to be able to visualize whatever the likelihood function represents independent of the probabilistic model.
+  - This would be useful to be able to visualize whatever the likelihood function represents independent of the probabilistic model.
+  - 
+### What would I want to have to do StonePark example as flowless as possible? 
+- A way to save the model already with the correct **nugget effects** instead of having to "load nuggets"
+- Easy pykeops setting 
+- Likelihood function plug and play
+- Some sort of pre-analysis if the model looks good. Maybe some metrics about how long it takes to run the model. 
+- Some sort of easy set-up for having a reduced number of priors

examples/examples/_aux_func_gravity_inversion.py

@@ -1,6 +1,13 @@
 import numpy as np
 from sklearn.linear_model import LinearRegression
 
+import numpy as np
+
+import gempy as gp
+import xarray as xr
+from vector_geology.model_building_functions import optimize_nuggets_for_group
+
+
 def calculate_scale_shift(a: np.ndarray, b: np.ndarray) -> tuple:
     # Reshape arrays for sklearn
     a_reshaped = a.reshape(-1, 1)
@@ -22,17 +29,10 @@ def gaussian_kernel(locations, length_scale, variance):
     locations = torch.tensor(locations.values)
     distance_squared = torch.cdist(locations, locations, p=2).pow(2)
     # Compute the covariance matrix using the Gaussian kernel
-    covariance_matrix = variance * torch.exp(-0.5 * distance_squared / length_scale**2)
+    covariance_matrix = variance * torch.exp(-0.5 * distance_squared / length_scale ** 2)
     return covariance_matrix
 
 
-import numpy as np
-
-import gempy as gp
-import xarray as xr
-from vector_geology.model_building_functions import optimize_nuggets_for_group
-
-
 def initialize_geo_model(structural_elements: list[gp.data.StructuralElement], extent: list[float],
                          topography: xr.Dataset, load_nuggets: bool = False
                          ) -> gp.data.GeoModel:

examples/examples/_aux_func_gravity_inversion_II.py

@@ -0,0 +1,94 @@
+import os
+import xarray as xr
+import numpy as np
+import gempy as gp
+import subsurface
+
+
+def process_file(filename, global_extent, color_generator: gp.data.ColorsGenerator):
+    dataset: xr.Dataset = xr.open_dataset(filename)
+    file_extent = _calculate_extent(dataset)
+    global_extent = _update_global_extent(global_extent, file_extent)
+
+    base, ext = os.path.splitext(filename)
+    base = os.path.basename(base)
+    structural_element = _extract_surface_points_and_orientations(dataset, base, color_generator)
+
+    return structural_element, global_extent
+
+
+# Function to calculate the extent (ROI) for a single dataset
+def _calculate_extent(dataset):
+    x_coord, y_coord, z_coord = dataset.vertex[:, 0], dataset.vertex[:, 1], dataset.vertex[:, 2]
+    return [x_coord.min().values, x_coord.max().values, y_coord.min().values, y_coord.max().values, z_coord.min().values, z_coord.max().values]
+
+
+# Function to extract surface points and orientations from the dataset
+
+
+def _extract_surface_points_and_orientations(dataset, name, color_generator) -> gp.data.StructuralElement:
+    name = name.replace("Stonepark_", "")
+
+    # Extract surface points and orientations
+    unstruct = subsurface.UnstructuredData(dataset)
+    ts = subsurface.TriSurf(mesh=unstruct)
+    triangulated_mesh = subsurface.visualization.to_pyvista_mesh(ts)
+    # print(name)
+    # subsurface.visualization.pv_plot([triangulated_mesh])
+
+    # Decimate the mesh and compute normals
+    decimated_mesh = triangulated_mesh.decimate_pro(0.95)
+    normals = decimated_mesh.compute_normals(point_normals=False, cell_normals=True, consistent_normals=True)
+    normals_array = np.array(normals.cell_data["Normals"])
+
+    # Filter points
+    vertex_sp = normals.cell_centers().points
+    vertex_grads = normals_array
+    positive_z = vertex_grads[:, 2] > 0
+    largest_component_z = np.all(vertex_grads[:, 2:3] >= np.abs(vertex_grads[:, :2]), axis=1)
+    filter_mask = np.logical_and(positive_z, largest_component_z)
+
+    surface_points_xyz = vertex_sp[filter_mask]
+    orientations_gxyz = vertex_grads[filter_mask]
+    nuggets = np.ones(len(surface_points_xyz)) * 0.000001
+
+    # Create SurfacePointsTable and OrientationsTable
+    every = 10
+    surface_points = gp.data.SurfacePointsTable.from_arrays(
+        x=surface_points_xyz[::every, 0],
+        y=surface_points_xyz[::every, 1],
+        z=surface_points_xyz[::every, 2],
+        names=name,
+        nugget=nuggets[::every]
+    )
+
+    every = 10
+    orientations = gp.data.OrientationsTable.from_arrays(
+        x=surface_points_xyz[::every, 0],
+        y=surface_points_xyz[::every, 1],
+        z=surface_points_xyz[::every, 2],
+        G_x=orientations_gxyz[::every, 0],
+        G_y=orientations_gxyz[::every, 1],
+        G_z=orientations_gxyz[::every, 2],
+        nugget=4,
+        names=name
+    )
+
+    structural_element = gp.data.StructuralElement(
+        name=name,
+        surface_points=surface_points,
+        orientations=orientations,
+        color=next(color_generator)
+    )
+
+    return structural_element
+
+
+# Function to update the global extent based on each file's extent
+def _update_global_extent(global_extent, file_extent):
+    if not global_extent:
+        return file_extent
+    else:
+        return [min(global_extent[0], file_extent[0]), max(global_extent[1], file_extent[1]),
+                min(global_extent[2], file_extent[2]), max(global_extent[3], file_extent[3]),
+                min(global_extent[4], file_extent[4]), max(global_extent[5], file_extent[5])]

examples/examples/gravity_inversion.py

@@ -1,4 +1,3 @@
-
 """
 Probabilistic Inversion Example: Geological Model
 --------------------------------------------------
@@ -22,12 +21,12 @@
 import gempy as gp
 import gempy_engine
 import gempy_viewer as gpv
+
+from examples.examples._aux_func_gravity_inversion import calculate_scale_shift, gaussian_kernel, initialize_geo_model, setup_geophysics 
+from examples.examples._aux_func_gravity_inversion_II import process_file
+
 from gempy_engine.core.backend_tensor import BackendTensor
 from gempy_probability.modules.plot.plot_posterior import default_red, default_blue
-from vector_geology.bayesian_helpers import calculate_scale_shift, gaussian_kernel
-from vector_geology.model_1_builder import initialize_geo_model, setup_geophysics
-from vector_geology.omf_to_gempy import process_file
-
 
 # %%
 # Config
@@ -65,7 +64,7 @@
     structural_elements=structural_elements,
     extent=(np.array(global_extent)),
     topography=(xr.open_dataset(os.path.join(path, "Topography.nc"))),
-    load_nuggets=True
+    load_nuggets=False
 )
 
 # %%
@@ -138,6 +137,7 @@
     densities=prior_tensor,
 )
 
+
 # %%
 # Define the Pyro probabilistic model for inversion
 def model(y_obs_list, interpolation_input):
@@ -168,6 +168,7 @@ def model(y_obs_list, interpolation_input):
         obs=y_obs_list
     )
 
+
 # %%
 # Prepare observed data for Pyro model and optimize mesh settings
 y_obs_list = torch.tensor(adapted_observed_grav.values).view(1, 17)
@@ -216,7 +217,7 @@ def model(y_obs_list, interpolation_input):
 )
 plt.show()
 
-#%%
+# %%
 az.plot_density(
     data=[data.posterior_predictive, data.prior_predictive],
     shade=.9,