.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_probabilistic/06_enmap_inversion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_probabilistic_06_enmap_inversion.py: Bayesian EnMap Inversion: Categorical Likelihood and Ordinal Probabilities ========================================================================== This tutorial demonstrates how to perform Bayesian inversion using surface lithological classifications (e.g., from EnMap satellite data) as observations. Unlike gravity inversion which uses continuous values, EnMap inversion deals with discrete categorical labels. **What are we doing?** We're inferring geological parameters (like layer orientations) by comparing the predicted lithology at the surface with observed lithological labels derived from remote sensing data. **Why Bayesian Inference?** Bayesian inference allows us to: 1. **Quantify Uncertainty**: Instead of a single "best-fit" model, we get a posterior distribution representing all models consistent with both our prior knowledge and the EnMap observations. 2. **Incorporate Prior Knowledge**: We can formally include geological information (e.g., typical layer dips) through prior distributions. 3. **Handle Classification Errors**: The probabilistic approach naturally accounts for potential misclassifications in the remote sensing data. **The Bayesian Framework** We seek the posterior distribution of geological parameters :math:`\theta` given the observed EnMap labels :math:`y`: .. math:: p(\theta|y) = \frac{p(y|\theta) p(\theta)}{p(y)} Where: - :math:`p(\theta)` is the **prior**, encoding our initial knowledge. - :math:`p(y|\theta)` is the **likelihood**, describing the probability of the labels given a specific geological configuration. - :math:`p(y)` is the **evidence** (marginal likelihood). **The Forward Model** The relationship between parameters and observations is given by: .. math:: y = f(\theta) + \epsilon Where :math:`f(\theta)` is the GemPy geological model that predicts rock types at the surface, and :math:`\epsilon` represents the observation noise or classification uncertainty. **The Likelihood: Categorical with Ordinal Probabilities** Since the observed data consists of discrete rock type labels (e.g., 0 for Unit A, 1 for Unit B), we use a **Categorical likelihood**. However, GemPy's forward model is differentiable and produces a continuous scalar field. To bridge this gap, we use an **Ordinal Probabilities** approach. 1. **Scalar Field to Probabilities**: We take the continuous scalar field value at the surface. 2. **Sigmoid Boundaries**: We define boundaries (thresholds) between units. We use sigmoid functions centered at these boundaries to compute the probability of a point belonging to each unit. .. math:: P(\text{below boundary}) = \sigma\left(\frac{\text{boundary} - \text{scalar}}{\text{temperature}}\right) 3. **Softmax-like behavior**: The difference between cumulative probabilities at successive boundaries gives the probability mass for each discrete class. 4. **Temperature**: A 'temperature' parameter controls the sharpness of these probabilistic transitions. Lower temperature means sharper, more deterministic boundaries. **Why this approach?** - **Differentiability**: The sigmoid-based mapping is fully differentiable, which is crucial for Gradient-based MCMC methods like **NUTS** (No-U-Turn Sampler). - **Uncertainty**: It naturally accounts for classification uncertainty near geological interfaces. .. GENERATED FROM PYTHON SOURCE LINES 83-88 Import Libraries ---------------- We load the required scientific packages, GemPy ecosystem, and Pyro for probabilistic programming. .. GENERATED FROM PYTHON SOURCE LINES 88-123 .. code-block:: Python import os import numpy as np import matplotlib.pyplot as plt import torch import pyro import pyro.distributions as dist import arviz as az import inspect from pathlib import Path # GemPy ecosystem import gempy as gp import gempy_probability as gpp import gempy_viewer as gpv # Configuration and helpers from mineye.config import paths from gempy_probability.core.samplers_data import NUTSConfig from gempy_probability.modules.plot.plot_posterior import default_red, default_blue from mineye.GeoModel.model_one.probabilistic_model import _modify_orientations from mineye.GeoModel.model_one.visualization import ( gempy_viz, plot_probability_heatmap, compute_probability_density_fields, plot_many_observed_vs_forward ) from mineye.GeoModel.model_one.probabilistic_model_likelihoods import enmap_likelihood_fn # Set random seeds for reproducibility seed = 4003 pyro.set_rng_seed(seed) torch.manual_seed(seed) np.random.seed(seed) .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH .. GENERATED FROM PYTHON SOURCE LINES 124-134 Step 1: Model Configuration --------------------------- We define the spatial extent of our study area and the resolution of the geological model. **Model Extent**: - X: -707521 to -675558 - Y: 4526832 to 4551949 - Z: -500 to 505 (meters) .. GENERATED FROM PYTHON SOURCE LINES 134-138 .. code-block:: Python # Define Model Extent (matching the gravity example area) extent = [-707521, -675558, 4526832, 4551949, -500, 505] .. GENERATED FROM PYTHON SOURCE LINES 139-149 Step 2: Load EnMap Observations ------------------------------- **What is EnMap data?** EnMap (Environmental Mapping and Analysis Program) provides hyperspectral satellite imagery. In this example, we use a lithological classification derived from such data. **Units**: Discrete categorical labels (0, 1, 2...). .. GENERATED FROM PYTHON SOURCE LINES 149-166 .. code-block:: Python # Define paths to find the central_xyz files in the same folder as this script current_dir = Path(inspect.getfile(inspect.currentframe())).parent.resolve() xyz_path = os.path.join(current_dir, 'central_xyz.npy') labels_path = os.path.join(current_dir, 'central_labels.npy') # Load EnMap data xyz_enmap = np.load(xyz_path) labels_enmap = np.load(labels_path) labels_enmap[labels_enmap == 2] = 1 # Normalize labels to binary for this example print(f"\nEnMap observations:") print(f" Number of measurements: {len(xyz_enmap)}") print(f" Classes present: {np.unique(labels_enmap)}") print(f" Class 0 count: {np.sum(labels_enmap == 0)}") print(f" Class 1 count: {np.sum(labels_enmap == 1)}") .. rst-class:: sphx-glr-script-out .. code-block:: none EnMap observations: Number of measurements: 57 Classes present: [0 1] Class 0 count: 27 Class 1 count: 30 .. GENERATED FROM PYTHON SOURCE LINES 167-171 Step 3: Initial Geological Model Setup -------------------------------------- We create the base GemPy model using geological orientations and points. .. GENERATED FROM PYTHON SOURCE LINES 171-196 .. code-block:: Python # Create Initial GemPy Model simple_geo_model = gp.create_geomodel( project_name='enmap_inversion', extent=extent, refinement=5, importer_helper=gp.data.ImporterHelper( path_to_orientations=paths.get_orientations_path(), path_to_surface_points=paths.get_points_path(), ) ) gp.map_stack_to_surfaces( gempy_model=simple_geo_model, mapping_object={"Tournaisian_Plutonites": ["Tournaisian Plutonites"]} ) # Compute initial model gp.compute_model(simple_geo_model) # Visualize initial model gpv.plot_2d(simple_geo_model, direction='z', show_data=True) plt.title("Initial Geological Model (Surface)") plt.show() .. image-sg:: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_001.png :alt: Initial Geological Model (Surface) :srcset: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces /home/leguark/.venv/2025/lib/python3.13/site-packages/torch/utils/_device.py:103: UserWarning: Using torch.cross without specifying the dim arg is deprecated. Please either pass the dim explicitly or simply use torch.linalg.cross. The default value of dim will change to agree with that of linalg.cross in a future release. (Triggered internally at /pytorch/aten/src/ATen/native/Cross.cpp:63.) return func(*args, **kwargs) .. GENERATED FROM PYTHON SOURCE LINES 197-202 Step 4: Configure Custom Grid for Observations ---------------------------------------------- We set a custom grid at the exact XYZ locations where we have EnMap labels. This ensures GemPy predicts rock types exactly where we have observations. .. GENERATED FROM PYTHON SOURCE LINES 202-210 .. code-block:: Python gp.set_custom_grid(simple_geo_model.grid, xyz_enmap) gp.set_active_grid( grid=simple_geo_model.grid, grid_type=[simple_geo_model.grid.GridTypes.CUSTOM], reset=True ) .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.OCTREE|CUSTOM|NONE Active grids: GridTypes.CUSTOM|NONE Grid(values=array([[-6.86311507e+05, 4.53422704e+06, 1.31536606e+02], [-6.84340294e+05, 4.53983741e+06, 3.46467194e+02], [-6.96432927e+05, 4.54468963e+06, 3.20259308e+02], [-6.91163723e+05, 4.54302168e+06, 2.81383545e+02], [-6.89533682e+05, 4.53085323e+06, 3.19174347e+02], [-6.95447321e+05, 4.53923089e+06, 4.63476868e+02], [-6.91504895e+05, 4.53479566e+06, 2.90548767e+02], [-6.81648830e+05, 4.53252118e+06, 1.20361305e+02], [-6.94916609e+05, 4.53570545e+06, 3.06075592e+02], [-7.05151754e+05, 4.54359030e+06, 3.37633240e+02], [-7.02460290e+05, 4.54700201e+06, 4.37251770e+02], [-7.00868156e+05, 4.54920067e+06, 5.18035645e+02], [-6.98025061e+05, 4.53892762e+06, 4.07593109e+02], [-7.01626315e+05, 4.54396938e+06, 3.42894135e+02], [-6.98138785e+05, 4.53225582e+06, 3.96869446e+02], [-7.01702131e+05, 4.53005716e+06, 3.33882172e+02], [-7.05341293e+05, 4.53399959e+06, 4.36949951e+02], [-7.04279871e+05, 4.53619825e+06, 4.22504150e+02], [-6.78957366e+05, 4.53585708e+06, 3.16025940e+02], [-6.77668496e+05, 4.53824528e+06, 3.09397583e+02], [-6.82558621e+05, 4.52861666e+06, 2.47029175e+02], [-6.78995274e+05, 4.54571314e+06, 3.88674500e+02], [-7.04962214e+05, 4.53877599e+06, 3.19339355e+02], [-6.80094604e+05, 4.52941273e+06, 2.73347809e+02], [-6.83657951e+05, 4.54529616e+06, 3.62819458e+02], [-6.96849914e+05, 4.52865457e+06, 2.98927002e+02], [-6.86576862e+05, 4.54685038e+06, 4.12333984e+02], [-7.01323052e+05, 4.53976160e+06, 3.07448608e+02], [-6.89950669e+05, 4.53972369e+06, 2.46116455e+02], [-6.80056696e+05, 4.54321122e+06, 3.16967072e+02], [-6.95030333e+05, 4.54817716e+06, 2.64292694e+02], [-7.03673344e+05, 4.54525825e+06, 4.04202576e+02], [-6.95636860e+05, 4.53134603e+06, 2.56243866e+02], [-7.05038030e+05, 4.54040603e+06, 3.10769226e+02], [-7.01019788e+05, 4.53657733e+06, 2.82767090e+02], [-6.79753433e+05, 4.53468193e+06, 2.85022522e+02], [-6.99996274e+05, 4.53172511e+06, 4.26705627e+02], [-6.97342718e+05, 4.53616034e+06, 2.89172455e+02], [-6.89761129e+05, 4.54499289e+06, 3.07890228e+02], [-6.79412261e+05, 4.53733549e+06, 3.21680176e+02], [-6.94006819e+05, 4.53445449e+06, 3.16392212e+02], [-6.84302386e+05, 4.54324913e+06, 2.76777466e+02], [-6.85325900e+05, 4.54582687e+06, 3.95587402e+02], [-6.86842218e+05, 4.54355239e+06, 3.20481750e+02], [-6.95902216e+05, 4.54203607e+06, 1.92909485e+02], [-7.03445896e+05, 4.53172511e+06, 3.84196045e+02], [-7.05379201e+05, 4.52876829e+06, 2.67404938e+02], [-6.94764978e+05, 4.52854085e+06, 2.83653137e+02], [-6.99693010e+05, 4.54552360e+06, 3.34157135e+02], [-6.99768826e+05, 4.54787390e+06, 2.46570984e+02], [-7.04393595e+05, 4.53471984e+06, 4.75726532e+02], [-6.99996274e+05, 4.54290796e+06, 3.11612671e+02], [-6.87676193e+05, 4.52945064e+06, 2.93266357e+02], [-6.84567741e+05, 4.52929901e+06, 3.24543060e+02], [-6.83278871e+05, 4.53161139e+06, 2.89459503e+02], [-7.05682465e+05, 4.54836670e+06, 3.66013245e+02], [-6.88282720e+05, 4.53183884e+06, 1.89883438e+02]]), length=array([], dtype=float64), _octree_grid=RegularGrid(resolution=array([512, 400, 32]), extent=array([-7.075210e+05, -6.755580e+05, 4.526832e+06, 4.551949e+06, -5.000000e+02, 5.050000e+02]), values=array([[-7.07489786e+05, 4.52686340e+06, -4.84296875e+02], [-7.07489786e+05, 4.52686340e+06, -4.52890625e+02], [-7.07489786e+05, 4.52686340e+06, -4.21484375e+02], ..., [-6.75589214e+05, 4.55191760e+06, 4.26484375e+02], [-6.75589214e+05, 4.55191760e+06, 4.57890625e+02], [-6.75589214e+05, 4.55191760e+06, 4.89296875e+02]], shape=(6553600, 3)), mask_topo=array([], shape=(0, 3), dtype=bool), _transform=None, _base_resolution=array([32, 25, 2])), _dense_grid=None, _custom_grid=CustomGrid(values=array([[-6.86311507e+05, 4.53422704e+06, 1.31536606e+02], [-6.84340294e+05, 4.53983741e+06, 3.46467194e+02], [-6.96432927e+05, 4.54468963e+06, 3.20259308e+02], [-6.91163723e+05, 4.54302168e+06, 2.81383545e+02], [-6.89533682e+05, 4.53085323e+06, 3.19174347e+02], [-6.95447321e+05, 4.53923089e+06, 4.63476868e+02], [-6.91504895e+05, 4.53479566e+06, 2.90548767e+02], [-6.81648830e+05, 4.53252118e+06, 1.20361305e+02], [-6.94916609e+05, 4.53570545e+06, 3.06075592e+02], [-7.05151754e+05, 4.54359030e+06, 3.37633240e+02], [-7.02460290e+05, 4.54700201e+06, 4.37251770e+02], [-7.00868156e+05, 4.54920067e+06, 5.18035645e+02], [-6.98025061e+05, 4.53892762e+06, 4.07593109e+02], [-7.01626315e+05, 4.54396938e+06, 3.42894135e+02], [-6.98138785e+05, 4.53225582e+06, 3.96869446e+02], [-7.01702131e+05, 4.53005716e+06, 3.33882172e+02], [-7.05341293e+05, 4.53399959e+06, 4.36949951e+02], [-7.04279871e+05, 4.53619825e+06, 4.22504150e+02], [-6.78957366e+05, 4.53585708e+06, 3.16025940e+02], [-6.77668496e+05, 4.53824528e+06, 3.09397583e+02], [-6.82558621e+05, 4.52861666e+06, 2.47029175e+02], [-6.78995274e+05, 4.54571314e+06, 3.88674500e+02], [-7.04962214e+05, 4.53877599e+06, 3.19339355e+02], [-6.80094604e+05, 4.52941273e+06, 2.73347809e+02], [-6.83657951e+05, 4.54529616e+06, 3.62819458e+02], [-6.96849914e+05, 4.52865457e+06, 2.98927002e+02], [-6.86576862e+05, 4.54685038e+06, 4.12333984e+02], [-7.01323052e+05, 4.53976160e+06, 3.07448608e+02], [-6.89950669e+05, 4.53972369e+06, 2.46116455e+02], [-6.80056696e+05, 4.54321122e+06, 3.16967072e+02], [-6.95030333e+05, 4.54817716e+06, 2.64292694e+02], [-7.03673344e+05, 4.54525825e+06, 4.04202576e+02], [-6.95636860e+05, 4.53134603e+06, 2.56243866e+02], [-7.05038030e+05, 4.54040603e+06, 3.10769226e+02], [-7.01019788e+05, 4.53657733e+06, 2.82767090e+02], [-6.79753433e+05, 4.53468193e+06, 2.85022522e+02], [-6.99996274e+05, 4.53172511e+06, 4.26705627e+02], [-6.97342718e+05, 4.53616034e+06, 2.89172455e+02], [-6.89761129e+05, 4.54499289e+06, 3.07890228e+02], [-6.79412261e+05, 4.53733549e+06, 3.21680176e+02], [-6.94006819e+05, 4.53445449e+06, 3.16392212e+02], [-6.84302386e+05, 4.54324913e+06, 2.76777466e+02], [-6.85325900e+05, 4.54582687e+06, 3.95587402e+02], [-6.86842218e+05, 4.54355239e+06, 3.20481750e+02], [-6.95902216e+05, 4.54203607e+06, 1.92909485e+02], [-7.03445896e+05, 4.53172511e+06, 3.84196045e+02], [-7.05379201e+05, 4.52876829e+06, 2.67404938e+02], [-6.94764978e+05, 4.52854085e+06, 2.83653137e+02], [-6.99693010e+05, 4.54552360e+06, 3.34157135e+02], [-6.99768826e+05, 4.54787390e+06, 2.46570984e+02], [-7.04393595e+05, 4.53471984e+06, 4.75726532e+02], [-6.99996274e+05, 4.54290796e+06, 3.11612671e+02], [-6.87676193e+05, 4.52945064e+06, 2.93266357e+02], [-6.84567741e+05, 4.52929901e+06, 3.24543060e+02], [-6.83278871e+05, 4.53161139e+06, 2.89459503e+02], [-7.05682465e+05, 4.54836670e+06, 3.66013245e+02], [-6.88282720e+05, 4.53183884e+06, 1.89883438e+02]])), _topography=None, _sections=None, _centered_grid=None, _transform=None, _octree_levels=-1) .. GENERATED FROM PYTHON SOURCE LINES 211-216 Step 5: Baseline Forward Model ------------------------------ We compute the model at the custom grid locations to establish a baseline comparison with the observations. .. GENERATED FROM PYTHON SOURCE LINES 216-227 .. code-block:: Python # Compute model on custom grid sol = gp.compute_model(simple_geo_model) # The results at custom grid are stored in solutions output_center = sol.octrees_output[0].last_output_center scalar_field_at_custom_grid = output_center.exported_fields.scalar_field[output_center.grid.custom_grid_slice] print(f"\nBaseline Forward Model:") print(f" Scalar field range: {scalar_field_at_custom_grid.min():.2f} to {scalar_field_at_custom_grid.max():.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Baseline Forward Model: Scalar field range: -0.04 to 0.05 .. GENERATED FROM PYTHON SOURCE LINES 228-247 Step 6: Define Prior Distributions ----------------------------------- **Prior Knowledge in Bayesian Inference** Priors encode our geological knowledge *before* seeing the data. Here we specify our uncertainty about the dip of the geological layers. **Prior on Dips** We use a Normal distribution: .. math:: \text{dips} \sim \mathcal{N}(10^\circ, 10^\circ) **Interpretation**: - **Mean (μ)**: 10 degrees. Our best guess for the regional dip. - **Std (σ)**: 10 degrees. Represents significant uncertainty in our knowledge. .. GENERATED FROM PYTHON SOURCE LINES 247-258 .. code-block:: Python model_priors = { 'dips': dist.Normal( loc=(torch.ones(simple_geo_model.orientations_copy.xyz.shape[0]) * 10), scale=torch.tensor(10, dtype=torch.float64) ).to_event(1) } def set_priors_enmap(samples, geo_model): return _modify_orientations(samples=samples, geo_model=geo_model, key="dips") .. GENERATED FROM PYTHON SOURCE LINES 259-267 Step 7: Inspecting the Likelihood Function ------------------------------------------ **Understanding the Likelihood Function** The likelihood function is crucial for connecting model predictions to observations. For EnMap data, we use an ordinal classification likelihood. Let's inspect what it does internally: .. GENERATED FROM PYTHON SOURCE LINES 267-271 .. code-block:: Python # Display the function signature and docstring help(enmap_likelihood_fn) .. rst-class:: sphx-glr-script-out .. code-block:: none Help on function enmap_likelihood_fn in module mineye.GeoModel.model_one.probabilistic_model_likelihoods: enmap_likelihood_fn(solutions: gempy_engine.core.data.solutions.Solutions) .. GENERATED FROM PYTHON SOURCE LINES 272-276 **Full Source Code** For complete transparency, here's the full implementation of the EnMap likelihood function: .. GENERATED FROM PYTHON SOURCE LINES 276-281 .. code-block:: Python print("Likelihood function source code:") print("=" * 60) print(inspect.getsource(enmap_likelihood_fn)) .. rst-class:: sphx-glr-script-out .. code-block:: none Likelihood function source code: ============================================================ def enmap_likelihood_fn(solutions: gp.data.Solutions): output_center: "gempy_engine.core.data.interp_output.InterpOutput" = solutions.octrees_output[0].last_output_center scalar_field_at_custom_grid = output_center.exported_fields.scalar_field[output_center.grid.custom_grid_slice] if not isinstance(scalar_field_at_custom_grid, torch.Tensor): scalar_field_at_custom_grid = torch.tensor(scalar_field_at_custom_grid, dtype=torch.float64) # DEFINE YOUR BOUNDARIES HERE # These should match the isovalues of your interfaces. # If your model is normalized, they might be roughly integers or fixed floats. # Example: If Unit 0 is < 1.0, Unit 1 is 1.0 to 2.0, Unit 2 is > 2.0 # boundaries = torch.tensor([1.0, 2.0], dtype=torch.float64) boundaries = torch.tensor(solutions.scalar_field_at_surface_points, dtype=torch.float64) # Calculate probabilities probs = _get_ordinal_probs(scalar_field_at_custom_grid, boundaries, temperature=0.1) # 2. Save the class probabilities (The "Confidence") # This will be shape (n_points, n_classes) pyro.deterministic("probs_pred", probs.detach()) # 3. (Optional) Save the boundaries if they change per iteration pyro.deterministic("boundaries_pred", boundaries) # Use Categorical directly with probabilities (no need for logits) return dist.Categorical(probs=probs).to_event(1) .. GENERATED FROM PYTHON SOURCE LINES 282-293 Step 8: Probabilistic Model Creation ------------------------------------ We combine the priors, the geological model update function, and the likelihood function into a single Pyro model. **Internal Workflow**: 1. Sample :math:`\theta` (dips) from priors. 2. Update GemPy model with sampled :math:`\theta`. 3. Compute forward model (predict scalar fields at surface). 4. Apply `enmap_likelihood_fn` to compute :math:`p(data|\theta)`. .. GENERATED FROM PYTHON SOURCE LINES 293-305 .. code-block:: Python # Create the Pyro model prob_model = gpp.make_gempy_pyro_model( priors=model_priors, set_interp_input_fn=set_priors_enmap, likelihood_fn=enmap_likelihood_fn, obs_name="EnMap Labels" ) # Prepare observed data as a torch tensor labels_enmap_tensor = torch.tensor(labels_enmap, dtype=torch.float64) .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH .. GENERATED FROM PYTHON SOURCE LINES 306-312 Step 9: Prior Predictive Check ------------------------------ Before running expensive inference, we check if our priors allow for the observed rock types. If the prior predictions don't even cover the observations, the model might be poorly specified. .. GENERATED FROM PYTHON SOURCE LINES 312-334 .. code-block:: Python print("Running prior predictive...") prior_data = gpp.run_predictive( prob_model=prob_model, geo_model=simple_geo_model, y_obs_list=None, n_samples=100 ) # Plotting the prior distribution of labels vs observed pred_labels = prior_data.prior['EnMap Labels'][0] plt.figure(figsize=(8, 4)) plt.hist([pred_labels.values.flatten(), labels_enmap], label=['Prior Predicted', 'Observed Ground Truth'], bins=np.arange(3) - 0.5, density=True, rwidth=0.8) plt.xticks([0, 1]) plt.xlabel("Rock Type Label") plt.ylabel("Density") plt.title("Prior Predictive Check") plt.legend() plt.show() .. image-sg:: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_002.png :alt: Prior Predictive Check :srcset: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Running prior predictive... Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces Using sequential processing for 1 surfaces .. GENERATED FROM PYTHON SOURCE LINES 335-341 Step 10: MCMC Inference using NUTS ---------------------------------- We use the No-U-Turn Sampler (NUTS) to sample from the posterior distribution. NUTS is an efficient Hamiltonian Monte Carlo algorithm that avoids the need for manual tuning of step sizes. .. GENERATED FROM PYTHON SOURCE LINES 341-376 .. code-block:: Python RUN_SIMULATION = False if RUN_SIMULATION: print("Running NUTS inference...") data = gpp.run_nuts_inference( prob_model=prob_model, geo_model=simple_geo_model, y_obs_list=labels_enmap_tensor, config=NUTSConfig( num_samples=200, warmup_steps=200, target_accept_prob=0.65, max_tree_depth=5 ), plot_trace=True, run_posterior_predictive=True ) data.extend(prior_data) data.to_netcdf("arviz_data_enmap.nc") else: # Get the directory of the current file current_dir = Path(inspect.getfile(inspect.currentframe())).parent.resolve() data_path = current_dir / "arviz_data_enmap_feb2026.nc" if not data_path.exists(): raise FileNotFoundError( f"Data file not found at {data_path}. " f"Please run the simulation first with RUN_SIMULATION=True" ) # Read the data file data = az.from_netcdf(str(data_path)) print(f"Loaded inference results from {data_path}") .. rst-class:: sphx-glr-script-out .. code-block:: none Loaded inference results from /home/leguark/PycharmProjects/Mineye-Terranigma/examples/02_probabilistic_modeling/arviz_data_enmap_feb2026.nc .. GENERATED FROM PYTHON SOURCE LINES 377-382 Step 11: Posterior Analysis --------------------------- We visualize how the data has informed our knowledge of the layer dips and assess the reduction in uncertainty. .. GENERATED FROM PYTHON SOURCE LINES 382-404 .. code-block:: Python # Posterior Density Plot az.plot_density( data=[data, data.prior], var_names=["dips"], data_labels=["Posterior", "Prior"], colors=[default_red, default_blue], shade=0.2 ) plt.title("Update of Dip Distributions") plt.show() # Uncertainty Reduction Calculation prior_std = data.prior["dips"].values.std() post_std = data.posterior["dips"].values.std() reduction = (1 - post_std / prior_std) * 100 print(f"\nUncertainty Reduction in Dips: {reduction:.1f}%") # Visualize the resulting geological model uncertainty print("\nVisualizing geological uncertainty...") gempy_viz(simple_geo_model, data) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_003.png :alt: dips 0, Update of Dip Distributions :srcset: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_004.png :alt: Cell Number: mid Direction: y :srcset: /examples_probabilistic/images/sphx_glr_06_enmap_inversion_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Uncertainty Reduction in Dips: 4.6% Visualizing geological uncertainty... Active grids: GridTypes.OCTREE|NONE Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces /home/leguark/PycharmProjects/gempy_viewer/gempy_viewer/modules/plot_2d/drawer_contours_2d.py:38: UserWarning: The following kwargs were not used by contour: 'linewidth', 'contour_colors' contour_set = ax.contour( Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces .. GENERATED FROM PYTHON SOURCE LINES 405-443 Step 12: Summary and Conclusions -------------------------------- **Workflow Recap** 1. **Prior Definition**: Normal distribution on layer dips. 2. **Likelihood Construction**: Ordinal classification using sigmoids. 3. **Model Integration**: Pyro/GemPy bridge. 4. **Inference**: NUTS sampling. 5. **Validation**: Posterior density and geological visualization. **Key Insights**: - Categorical data can be used for inversion by mapping continuous fields to probabilities via sigmoid functions. - Differentiability is preserved, allowing for efficient gradient-based sampling. - Remote sensing observations can significantly reduce uncertainty in subsurface geometrical parameters. Final workflow visualization (ASCII art) ┌─────────────────────────────────────────────────────────────────┐ │ 1. PRIOR DEFINITION │ │ dips ~ Normal(10, 10) │ └─────────────────────────────────────────────────────────────────┘ ↓ ┌─────────────────────────────────────────────────────────────────┐ │ 2. PROBABILISTIC MODEL CONSTRUCTION │ │ Bridge GemPy and Pyro with enmap_likelihood_fn │ └─────────────────────────────────────────────────────────────────┘ ↓ ┌─────────────────────────────────────────────────────────────────┐ │ 3. MCMC INFERENCE (NUTS) │ │ Sample posterior p(dips | EnMap Labels) │ └─────────────────────────────────────────────────────────────────┘ ↓ ┌─────────────────────────────────────────────────────────────────┐ │ 4. POSTERIOR ANALYSIS │ │ Compare prior vs posterior distributions │ └─────────────────────────────────────────────────────────────────┘ sphinx_gallery_thumbnail_number = 2 .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 18.013 seconds) .. _sphx_glr_download_examples_probabilistic_06_enmap_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06_enmap_inversion.ipynb <06_enmap_inversion.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06_enmap_inversion.py <06_enmap_inversion.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06_enmap_inversion.zip <06_enmap_inversion.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_