.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples_probabilistic/02_error_propagation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_probabilistic_02_error_propagation.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_probabilistic_02_error_propagation.py:


Error Propagation in Geological Models
=======================================

This example demonstrates uncertainty quantification through error propagation.
We perturb a surface point's Z-coordinate and propagate the uncertainty through the model.

.. GENERATED FROM PYTHON SOURCE LINES 10-12

Import Libraries
----------------

.. GENERATED FROM PYTHON SOURCE LINES 12-34

.. code-block:: Python


    import numpy as np
    import gempy as gp
    import gempy_viewer as gpv
    import gempy_probability as gpp

    import torch
    import pyro
    import pyro.distributions as dist
    from pyro.distributions import Distribution

    import arviz as az
    from gempy_engine.core.backend_tensor import BackendTensor
    from gempy_engine.core.data.interpolation_input import InterpolationInput
    from gempy_probability.modules.plot.plot_gempy import plot_gempy

    # Set random seeds for reproducibility
    seed = 4003
    pyro.set_rng_seed(seed)
    torch.manual_seed(seed)
    np.random.seed(1234)








.. GENERATED FROM PYTHON SOURCE LINES 35-36

Import paths configuration

.. GENERATED FROM PYTHON SOURCE LINES 36-38

.. code-block:: Python

    from mineye.config import paths








.. GENERATED FROM PYTHON SOURCE LINES 39-41

Define Model Extent and Resolution
-----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 41-54

.. code-block:: Python


    min_x = -707521
    max_x = -675558
    min_y = 4526832
    max_y = 4551949
    max_z = 505
    model_depth = -500
    extent = [min_x, max_x, min_y, max_y, model_depth, max_z]

    # Model resolution: use octree with refinement level 5
    resolution = None
    refinement = 5








.. GENERATED FROM PYTHON SOURCE LINES 55-57

Get Data Paths
--------------

.. GENERATED FROM PYTHON SOURCE LINES 57-64

.. code-block:: Python


    mod_or_path = paths.get_orientations_path()
    mod_pts_path = paths.get_points_path()

    print(f"Orientations: {mod_or_path}")
    print(f"Points: {mod_pts_path}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Orientations: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/orientations_mod.csv
    Points: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/points_mod.csv




.. GENERATED FROM PYTHON SOURCE LINES 65-67

Create GemPy Geological Model
------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 67-79

.. code-block:: Python


    geo_model = gp.create_geomodel(
        project_name='error_propagation_model',
        extent=extent,
        refinement=refinement,
        resolution=resolution,
        importer_helper=gp.data.ImporterHelper(
            path_to_orientations=mod_or_path,
            path_to_surface_points=mod_pts_path,
        )
    )








.. GENERATED FROM PYTHON SOURCE LINES 80-82

Map Geological Units
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 82-90

.. code-block:: Python


    gp.map_stack_to_surfaces(
        gempy_model=geo_model,
        mapping_object={
            "Tournaisian_Plutonites": ["Tournaisian Plutonites"],
        }
    )






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">

            <table>
              <tr><td>Structural Groups:</td><td>
        <table style="border-left:1.2px solid black;>
          <tr><th colspan="2"><b>StructuralGroup:</b></th></tr>
          <tr><td>Name:</td><td>Tournaisian_Plutonites</td></tr>
          <tr><td>Structural Relation:</td><td>StackRelationType.ERODE</td></tr>
          <tr><td>Elements:</td><td>
        <table width="50%" style="border-left:15px solid #015482;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>Tournaisian Plutonites</td></tr>
        </table>
            </td></tr>
        </table>
            </td></tr>
              <tr><td>Fault Relations:</td><td><table style="border-collapse: collapse; table-layout: fixed;"><th></th><th style="transform: rotate(-35deg); height:150px; vertical-align: bottom; text-align: center;">Tournaisia...</th><tr><th>Tournaisian_Plutonites</th><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td></tr></table></td></tr>
              <tr><td></td><td>
            <table>
              <tr>
                <td><div style="display: inline-block; background-color: #527682; width: 20px; height: 20px; border: 1px solid black;"></div> True</td>
                <td><div style="display: inline-block; background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></div> False</td>
              </tr>
            </table>
            </td></tr>
            </table>
        
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 91-94

Switch to PyTorch Backend
--------------------------
Required for probabilistic modeling

.. GENERATED FROM PYTHON SOURCE LINES 94-98

.. code-block:: Python


    BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH)
    print("✓ Switched to PyTorch backend")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Setting Backend To: AvailableBackends.PYTORCH
    ✓ Switched to PyTorch backend




.. GENERATED FROM PYTHON SOURCE LINES 99-117

Define Prior Distribution
--------------------------

**Prior Predictive Sampling vs. Posterior Sampling**

* **Prior Predictive Sampling**: Generating data from the model using parameters drawn
  from the prior distribution. It answers: "What kind of models/data do our initial
  beliefs produce?" It's a way to check if the model and priors are reasonable before
  considering observations.

* **Posterior Sampling**: Generating parameters that are consistent with both our
  initial beliefs (priors) AND the observed data. It answers: "What parameters are
  most likely given the data we've seen?"

In this example, we focus on **Prior Predictive Sampling** to see how uncertainty
in input coordinates propagates to the final structural geometry.

We'll add uncertainty to the Z-coordinate of the first surface point.

.. GENERATED FROM PYTHON SOURCE LINES 117-139

.. code-block:: Python


    from mineye.config.example_parameters import TharsisModelConfig
    FORMATION_COLORS = TharsisModelConfig.TharsisDataProcessingConfig.FORMATION_COLORS

    def modify_z_for_surface_point1(
        samples: dict[str, Distribution],
        geo_model: gp.data.GeoModel,
    ) -> InterpolationInput:
        """Modify the Z-coordinate of the first surface point based on samples."""
        prior_key = r'$\mu_{top}$'

        from gempy.modules.data_manipulation import interpolation_input_from_structural_frame
        interp_input = interpolation_input_from_structural_frame(geo_model)

        new_tensor: torch.Tensor = torch.index_put(
            input=interp_input.surface_points.sp_coords,
            indices=(torch.tensor([0]), torch.tensor([2])),
            values=(samples[prior_key])
        )
        interp_input.surface_points.sp_coords = new_tensor
        return interp_input








.. GENERATED FROM PYTHON SOURCE LINES 140-143

Set Prior Parameters
--------------------
Normal distribution around the original point location

.. GENERATED FROM PYTHON SOURCE LINES 143-156

.. code-block:: Python


    original_z = geo_model.surface_points_copy_transformed.xyz[0, 2]
    print(f"Original Z-coordinate: {original_z:.2f}")

    model_priors = {
        r'$\mu_{top}$': dist.Normal(
            loc=original_z,
            scale=torch.tensor(0.002, dtype=torch.float64)
        )
    }

    print(f"Prior: Normal(loc={original_z:.2f}, scale=0.001)")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Original Z-coordinate: 0.00
    Prior: Normal(loc=0.00, scale=0.001)




.. GENERATED FROM PYTHON SOURCE LINES 157-159

Create Probabilistic Model
---------------------------

.. GENERATED FROM PYTHON SOURCE LINES 159-169

.. code-block:: Python


    prob_model: gpp.GemPyPyroModel = gpp.make_gempy_pyro_model(
        priors=model_priors,
        set_interp_input_fn=modify_z_for_surface_point1,
        likelihood_fn=None,
        obs_name=None
    )

    print("✓ Probabilistic model created")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Setting Backend To: AvailableBackends.PYTORCH
    ✓ Probabilistic model created




.. GENERATED FROM PYTHON SOURCE LINES 170-173

Run Prior Predictive Sampling
------------------------------
Sample from the prior distribution and propagate through the model

.. GENERATED FROM PYTHON SOURCE LINES 173-187

.. code-block:: Python


    n_samples = 50
    print(f"Running {n_samples} prior predictive samples...")

    prior_inference_data: az.InferenceData = gpp.run_predictive(
        prob_model=prob_model,
        geo_model=geo_model,
        y_obs_list=[],
        n_samples=n_samples,
        plot_trace=True
    )

    print("✓ Prior predictive sampling complete")




.. image-sg:: /examples_probabilistic/images/sphx_glr_02_error_propagation_001.png
   :alt: $\mu_{top}$, $\mu_{top}$
   :srcset: /examples_probabilistic/images/sphx_glr_02_error_propagation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Running 50 prior predictive samples...
    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
    ✓ Prior predictive sampling complete




.. GENERATED FROM PYTHON SOURCE LINES 188-191

Visualize Uncertainty
---------------------
Plot the model with sampled surface points

.. GENERATED FROM PYTHON SOURCE LINES 191-228

.. code-block:: Python


    def update_model_for_plotting(geo_model: gp.data.GeoModel, sample_value: float, sample_idx: int):
        """Update model with a sampled Z-coordinate value."""
        xyz = np.zeros((1, 3))
        xyz[0, 2] = sample_value
        world_coord = geo_model.input_transform.apply_inverse(xyz)

        # Modify the surface point
        gp.modify_surface_points(
            geo_model=geo_model,
            slice=0,
            Z=world_coord[0, 2]
        )

    # Create base plot
    p2d = gpv.plot_2d(
        model=geo_model,
        show_topography=False,
        legend=False,
        show_lith=False,
        show_data=False,
        show=False,
        ve=5
    )

    # Overlay sampled models
    plot_gempy(
        geo_model=geo_model,
        n_samples=50,
        samples=(prior_inference_data.prior[r'$\mu_{top}$'].values[0, :]),
        update_model_fn=update_model_for_plotting,
        gempy_plot=p2d,
        contour_colors=[FORMATION_COLORS['Tournaisian Plutonites']]
    )

    print("✓ Visualization complete")




.. image-sg:: /examples_probabilistic/images/sphx_glr_02_error_propagation_002.png
   :alt: Cell Number: mid Direction: y
   :srcset: /examples_probabilistic/images/sphx_glr_02_error_propagation_002.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/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
    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
    ✓ Visualization complete




.. GENERATED FROM PYTHON SOURCE LINES 229-231

Summary Statistics
------------------

.. GENERATED FROM PYTHON SOURCE LINES 231-240

.. code-block:: Python


    samples = prior_inference_data.prior[r'$\mu_{top}$'].values[0, :]
    print(f"\nSampled Z-coordinate statistics:")
    print(f"  Mean: {samples.mean():.4f}")
    print(f"  Std: {samples.std():.4f}")
    print(f"  Min: {samples.min():.4f}")
    print(f"  Max: {samples.max():.4f}")
    print(f"  Original: {original_z:.4f}")

    # sphinx_gallery_thumbnail_number = 2



.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Sampled Z-coordinate statistics:
      Mean: 0.0010
      Std: 0.0017
      Min: -0.0023
      Max: 0.0046
      Original: 0.0004





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (1 minutes 0.859 seconds)


.. _sphx_glr_download_examples_probabilistic_02_error_propagation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 02_error_propagation.ipynb <02_error_propagation.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 02_error_propagation.py <02_error_propagation.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 02_error_propagation.zip <02_error_propagation.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_