.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "Tutorials\02_POC_create-MC-ensemble.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

        Click :ref:`here <sphx_glr_download_Tutorials_02_POC_create-MC-ensemble.py>`
        to download the full example code

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

.. _sphx_glr_Tutorials_02_POC_create-MC-ensemble.py:


02 - Monte Carlo simulation  
======================

The following tutorial will lead you through an example workflow on how to create a Monte Carlo simulation of 
geological models, meaning we will produce different geological geometries and also simulate their gravity
response.

.. GENERATED FROM PYTHON SOURCE LINES 10-13

Importing libraries
-------------------
First things first: let's import necessary libraries.

.. GENERATED FROM PYTHON SOURCE LINES 13-35

.. code-block:: default

    import warnings
    warnings.filterwarnings("ignore")

    # Importing GemPy
    import gempy as gp
    from gempy.assets import topology as tp
    from gempy.bayesian.fields import compute_prob, calculate_ie_masked
    from gempy.assets.geophysics import GravityPreprocessing

    # Importing auxilary libraries
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    plt.style.use('seaborn-talk')

    #import sys
    #sys.path.append('../../OpenWF/')
    #from aux_functions import log_progress

    # Check gempy version used for running the code
    print(f"Code run with GemPy version: {gp.__version__}")





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

 Out:

 .. code-block:: none

    Code run with GemPy version: 2.2.9




.. GENERATED FROM PYTHON SOURCE LINES 36-41

Model Initialization
--------------------

First, we import the base Proof-of-Concept model (POC-model from here on), which was generated in the previous example. Using the loading method of GemPy `gp.load_model()` directly loads the model's input, already set with fault relations, surfaces assigned to a stack (series), etc.
Only thing left is to recompile and run the model.

.. GENERATED FROM PYTHON SOURCE LINES 41-50

.. code-block:: default


    model_path = '../../models/2021-06-04_POC_base_model'

    geo_model = gp.load_model('POC_PCT_model', path=model_path,
                             recompile=False)

    # import DTM
    dtm = np.load('../../models/Graben_base_model_topography.npy')





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

 Out:

 .. code-block:: none

    Active grids: ['regular']
    Active grids: ['regular' 'topography']




.. GENERATED FROM PYTHON SOURCE LINES 51-52

Using the method `.get_additional_data()`, we can display a summary of model information and parameters, such as the kriging parameters.

.. GENERATED FROM PYTHON SOURCE LINES 52-55

.. code-block:: default


    geo_model.get_additional_data()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
          <th>values</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="9" valign="top">Structure</th>
          <th>isLith</th>
          <td>True</td>
        </tr>
        <tr>
          <th>isFault</th>
          <td>True</td>
        </tr>
        <tr>
          <th>number faults</th>
          <td>5</td>
        </tr>
        <tr>
          <th>number surfaces</th>
          <td>11</td>
        </tr>
        <tr>
          <th>number series</th>
          <td>10</td>
        </tr>
        <tr>
          <th>number surfaces per series</th>
          <td>[1, 1, 1, 1, 1, 3, 1, 1, 1, 0]</td>
        </tr>
        <tr>
          <th>len surfaces surface_points</th>
          <td>[6, 8, 12, 8, 6, 8, 14, 14, 26, 24, 24]</td>
        </tr>
        <tr>
          <th>len series surface_points</th>
          <td>[6, 8, 12, 8, 6, 36, 26, 24, 24, 0]</td>
        </tr>
        <tr>
          <th>len series orientations</th>
          <td>[2, 2, 6, 4, 2, 12, 16, 12, 10, 0]</td>
        </tr>
        <tr>
          <th rowspan="5" valign="top">Options</th>
          <th>dtype</th>
          <td>float64</td>
        </tr>
        <tr>
          <th>output</th>
          <td>geology</td>
        </tr>
        <tr>
          <th>theano_optimizer</th>
          <td>fast_compile</td>
        </tr>
        <tr>
          <th>device</th>
          <td>cpu</td>
        </tr>
        <tr>
          <th>verbosity</th>
          <td>None</td>
        </tr>
        <tr>
          <th rowspan="3" valign="top">Kriging</th>
          <th>range</th>
          <td>32190.8</td>
        </tr>
        <tr>
          <th>$C_o$</th>
          <td>2.46726e+07</td>
        </tr>
        <tr>
          <th>drift equations</th>
          <td>[3, 3, 3, 3, 3, 3, 3, 3, 3, 3]</td>
        </tr>
        <tr>
          <th rowspan="2" valign="top">Rescaling</th>
          <th>rescaling factor</th>
          <td>56916.7</td>
        </tr>
        <tr>
          <th>centers</th>
          <td>[14239.166495, 6900.0, -2581.92853]</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 56-57

Changing the kriging parameters affects the resulting models, e.g. the range represents the maximum correlation distance, or reducing the coefficient of correlation will yield a smoother, less "bumpy" model. For the POC-model, we set the `range` to 20000 and the correlation coefficient $C_o$ to 200000. Then we set up the interpolator, i.e. compile the functions which will calculate the scalar fields of our model surfaces.

.. GENERATED FROM PYTHON SOURCE LINES 57-75

.. code-block:: default


    # adapt kriging to the parameters of previous example
    # decrease the kriging range
    geo_model.modify_kriging_parameters('range', 20000.)
    geo_model.modify_kriging_parameters('$C_o$', 2e5)


    # Set the interpolator function
    # Create the theano model
    gp.set_interpolator(geo_model,
                             compile_theano=True,
                             theano_optimizer='fast_compile',
                             verbose=[],
                             update_kriging=False);

    # compute the model
    sol = gp.compute_model(geo_model, compute_mesh=True)





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

 Out:

 .. code-block:: none

    Compiling theano function...
    Level of Optimization:  fast_compile
    Device:  cpu
    Precision:  float64
    Number of faults:  5
    Compilation Done!
    Kriging values: 
                                              values
    range                                     20000
    $C_o$                                    200000
    drift equations  [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]




.. GENERATED FROM PYTHON SOURCE LINES 76-77

Now that the model is computed, lets have a look at a cross-section along the y-axis, so across the graben system:

.. GENERATED FROM PYTHON SOURCE LINES 77-81

.. code-block:: default


    gp.plot_2d(geo_model, cell_number=25, direction='y', show_data=False, show_topography=False,
              show_lith=True, show_results=True, show_boundaries=False);




.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_001.png
    :alt: Cell Number: 25 Direction: y
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none


    <gempy.plot.visualization_2d.Plot2D object at 0x0000017585940A60>



.. GENERATED FROM PYTHON SOURCE LINES 82-87

The two distinct domains in this model are directly visible: (i) the old graben system (extensional regime), covered by the (ii) thrusted, younger units.

Add Gravity grid
----------------
In the previous example, next to creating the model, we chose quasi-random locations for 15 gravity stations. The gravity signal of the base POC-model is simulated at these 15 stations. In the following workflows, we assume that these 15 stations were measured. So they serve as observed data for conditioning the MonteCarlo Ensemble of different geological geometries.

.. GENERATED FROM PYTHON SOURCE LINES 87-103

.. code-block:: default


    # In[7]:


    grav_stations = pd.read_csv('../../data/Data_for_MC/20210322_forw_grav_seed58.csv')
    station_coordinates = np.stack((grav_stations.X.values, 
                                    grav_stations.Y.values, 
                                    grav_stations.Z.values), axis=1)

    fig = plt.figure(figsize=[11,5])
    cb = plt.scatter(grav_stations['X'], grav_stations['Y'], c=grav_stations['grav'], 
               marker='s', s=90, cmap='viridis')
    plt.colorbar(cb, label='gravity')
    plt.ylabel('y [m]')
    plt.xlabel('x [m]');




.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_002.png
    :alt: 02 POC create MC ensemble
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none


    Text(0.5, 13.944444444444438, 'x [m]')



.. GENERATED FROM PYTHON SOURCE LINES 104-105

These stations are used for creating a centered grid around each station. The centered grid has an extent of 10 cells in x- and y-direction, and 15 cells in the z-direction.

.. GENERATED FROM PYTHON SOURCE LINES 105-112

.. code-block:: default




    geo_model.set_centered_grid(station_coordinates,  resolution = [10, 10, 15], radius=6000)
    g = GravityPreprocessing(geo_model.grid.centered_grid)
    tz = g.set_tz_kernel()





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

 Out:

 .. code-block:: none

    Active grids: ['regular' 'topography' 'centered']




.. GENERATED FROM PYTHON SOURCE LINES 113-118

We see that there are three active grids. On each, the gravity signal will be calculated. Of course, we can let it be calculated on each grid, but we may not need the information on e.g. the topography grid (which would for instance yield the geological map). 
So we can set only the centered grid to active, which speeds up the simulation.

**Note** that you'll need to model also the regular grid, if you plan to export the `lith_block` geological voxel model later on! 
As we want to also have the geometric changes in the lithological grid, we set `reset=False`. If we were to set it to `True`, only the 'centered' grid would be active.

.. GENERATED FROM PYTHON SOURCE LINES 118-121

.. code-block:: default


    geo_model.set_active_grid('centered', reset=False)





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

 Out:

 .. code-block:: none

    Active grids: ['regular' 'topography' 'centered']

    Grid Object. Values: 
    array([[  140.        ,   140.        , -6437.5       ],
           [  140.        ,   140.        , -6312.5       ],
           [  140.        ,   140.        , -6187.5       ],
           ...,
           [30888.89      ,  6285.71      , -3495.98176905],
           [30888.89      ,  6285.71      , -4948.49684561],
           [30888.89      ,  6285.71      , -6966.76      ]])



.. GENERATED FROM PYTHON SOURCE LINES 122-125

The centered grid will now be the only one where the model information is stored, meaning less computational time. Let's have a look how this comes in handy, when we start to modify the depth of units and calculate the gravity.

Before running the simulations, we need to assign densities to the rock units, otherwise it will raise an error.

.. GENERATED FROM PYTHON SOURCE LINES 125-131

.. code-block:: default


    # add densities - from abdelfettah 2014 and SAPHYR
    densities = [0, 0, 0, 0, 0, 2.466, 2.61, 2.53, 
                 2.61, 2.47, 2.55, 2.67]
    geo_model.add_surface_values(densities, ['density'])






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <style  type="text/css" >
        #T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col3 {
                background-color:  #5DA629;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col3 {
                background-color:  #5DA629;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col3 {
                background-color:  #015482;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col3 {
                background-color:  #015482;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col3 {
                background-color:  #015482;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col3 {
                background-color:  #dbdbac;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col3 {
                background-color:  #e588f3;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col3 {
                background-color:  #ff792b;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col3 {
                background-color:  #725c9a;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col3 {
                background-color:  #cfc199;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col3 {
                background-color:  #a5d490;
            }    #T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col3 {
                background-color:  #c7848f;
            }</style><table id="T_2f28591f_373c_11ec_9a60_00e04c6800ca" ><thead>    <tr>        <th class="blank level0" ></th>        <th class="col_heading level0 col0" >surface</th>        <th class="col_heading level0 col1" >series</th>        <th class="col_heading level0 col2" >order_surfaces</th>        <th class="col_heading level0 col3" >color</th>        <th class="col_heading level0 col4" >id</th>        <th class="col_heading level0 col5" >density</th>    </tr></thead><tbody>
                    <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row0" class="row_heading level0 row0" >9</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col0" class="data row0 col0" >Thrust1_south</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col1" class="data row0 col1" >Thrust1_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col2" class="data row0 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col3" class="data row0 col3" >#5DA629</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col4" class="data row0 col4" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow0_col5" class="data row0 col5" >0.000000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row1" class="row_heading level0 row1" >10</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col0" class="data row1 col0" >Thrust2_south</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col1" class="data row1 col1" >Thrust2_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col2" class="data row1 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col3" class="data row1 col3" >#5DA629</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col4" class="data row1 col4" >2</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow1_col5" class="data row1 col5" >0.000000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row2" class="row_heading level0 row2" >0</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col0" class="data row2 col0" >Fault2</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col1" class="data row2 col1" >Fault2_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col2" class="data row2 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col3" class="data row2 col3" >#015482</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col4" class="data row2 col4" >3</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow2_col5" class="data row2 col5" >0.000000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row3" class="row_heading level0 row3" >1</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col0" class="data row3 col0" >Fault5</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col1" class="data row3 col1" >Fault5_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col2" class="data row3 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col3" class="data row3 col3" >#015482</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col4" class="data row3 col4" >4</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow3_col5" class="data row3 col5" >0.000000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row4" class="row_heading level0 row4" >2</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col0" class="data row4 col0" >Fault6</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col1" class="data row4 col1" >Fault6_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col2" class="data row4 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col3" class="data row4 col3" >#015482</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col4" class="data row4 col4" >5</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow4_col5" class="data row4 col5" >0.000000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row5" class="row_heading level0 row5" >6</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col0" class="data row5 col0" >Tertiary</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col1" class="data row5 col1" >Post_tectonic_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col2" class="data row5 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col3" class="data row5 col3" >#dbdbac</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col4" class="data row5 col4" >6</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow5_col5" class="data row5 col5" >2.466000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row6" class="row_heading level0 row6" >8</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col0" class="data row6 col0" >Pink</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col1" class="data row6 col1" >Post_tectonic_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col2" class="data row6 col2" >2</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col3" class="data row6 col3" >#e588f3</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col4" class="data row6 col4" >7</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow6_col5" class="data row6 col5" >2.610000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row7" class="row_heading level0 row7" >7</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col0" class="data row7 col0" >Orange</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col1" class="data row7 col1" >Post_tectonic_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col2" class="data row7 col2" >3</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col3" class="data row7 col3" >#ff792b</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col4" class="data row7 col4" >8</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow7_col5" class="data row7 col5" >2.530000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row8" class="row_heading level0 row8" >5</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col0" class="data row8 col0" >Unconformity</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col1" class="data row8 col1" >Detachement</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col2" class="data row8 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col3" class="data row8 col3" >#725c9a</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col4" class="data row8 col4" >9</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow8_col5" class="data row8 col5" >2.610000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row9" class="row_heading level0 row9" >4</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col0" class="data row9 col0" >Upper-filling</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col1" class="data row9 col1" >Syn_tectonic_series2</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col2" class="data row9 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col3" class="data row9 col3" >#cfc199</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col4" class="data row9 col4" >10</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow9_col5" class="data row9 col5" >2.470000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row10" class="row_heading level0 row10" >3</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col0" class="data row10 col0" >Lower-filling</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col1" class="data row10 col1" >Pre_tectonic_series</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col2" class="data row10 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col3" class="data row10 col3" >#a5d490</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col4" class="data row10 col4" >11</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow10_col5" class="data row10 col5" >2.550000</td>
                </tr>
                <tr>
                            <th id="T_2f28591f_373c_11ec_9a60_00e04c6800calevel0_row11" class="row_heading level0 row11" >11</th>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col0" class="data row11 col0" >basement</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col1" class="data row11 col1" >Basement</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col2" class="data row11 col2" >1</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col3" class="data row11 col3" >#c7848f</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col4" class="data row11 col4" >12</td>
                            <td id="T_2f28591f_373c_11ec_9a60_00e04c6800carow11_col5" class="data row11 col5" >2.670000</td>
                </tr>
        </tbody></table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 132-135

MC Variation
------------
For varying the depth of units, we extract the indices of the units whose input points we want to modify. To guarantee that we always vary the original depth in each realization (and not the depth used in the previous realization), we first generate an initial-depth array, containing the original depth information of all input points:

.. GENERATED FROM PYTHON SOURCE LINES 135-138

.. code-block:: default


    Z_init = geo_model.surface_points.df['Z'].copy()








.. GENERATED FROM PYTHON SOURCE LINES 139-140

Having all the undisturbed depth values, we extract all surface points belonging to the units whose inputs we want to vary:

.. GENERATED FROM PYTHON SOURCE LINES 140-146

.. code-block:: default



    graben_lower = geo_model.surface_points.df.query("surface=='Lower-filling'")
    graben_middle = geo_model.surface_points.df.query("surface=='Upper-filling'")
    unconformity = geo_model.surface_points.df.query("surface=='Unconformity'")








.. GENERATED FROM PYTHON SOURCE LINES 147-148

Before running the Monte Carlo simulations, we set up the interpolator for a "fast-run", i.e. it optimizes runtime on cost of compilation time:

.. GENERATED FROM PYTHON SOURCE LINES 148-152

.. code-block:: default

    gp.set_interpolator(geo_model, output=['gravity'], 
                        theano_optimizer='fast_run', 
                        update_kriging=True)





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

 Out:

 .. code-block:: none

    Setting kriging parameters to their default values.
    Compiling theano function...
    Level of Optimization:  fast_run
    Device:  cpu
    Precision:  float64
    Number of faults:  5
    Compilation Done!
    Kriging values: 
                                              values
    range                                   32190.8
    $C_o$                               2.46726e+07
    drift equations  [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]

    <gempy.core.interpolator.InterpolatorModel object at 0x0000017586B9C340>



.. GENERATED FROM PYTHON SOURCE LINES 153-156

Now we are good to go and run the Monte Carlo simulation. In the following, we fix a numpy random number seed so that this MC-simulation is reproducible
Then, we create empty arrays and dictionaries for the lithologies and gravity, respectively. In a `for` loop, we then vary depths of interface points and
compute a model.

.. GENERATED FROM PYTHON SOURCE LINES 156-188

.. code-block:: default


    np.random.seed(1)
    # allocate array for lithology blocks
    lith_blocks = np.array([])
    # create a dictionary to store gravity of simulations
    grav = dict() 
    # get indices where the variable input points are
    Lgraben = list(graben_lower.index)
    Ugraben = list(graben_middle.index)
    Uncon = list(unconformity.index)
    Cindices = Lgraben + Ugraben + Uncon

    # set number of realizations
    n_iterations = 10

    for i in range(n_iterations):
        # vary surface points   
        Z_var = np.random.normal(0, 300, size=3)    
        Z_loc = np.hstack([Z_init[Lgraben] + Z_var[0],
                           Z_init[Ugraben] + Z_var[1],
                           Z_init[Uncon] + Z_var[2]])
        # apply variation to model
        geo_model.modify_surface_points(Cindices, Z=Z_loc)
        # re-compute model
        gp.compute_model(geo_model)
        # store lithologies ONLY THERE IF REGULAR GRID IS ACTIVE
        lith_blocks = np.append(lith_blocks, geo_model.solutions.lith_block)
        # store gravity
        grav[f"Real_{i}"] = geo_model.solutions.fw_gravity

    lith_blocks = lith_blocks.reshape(n_iterations, -1)








.. GENERATED FROM PYTHON SOURCE LINES 189-192

Export models and gravity
-------------------------
For post-processing of use in different software (e.g. numerical simulators for heat- and mass-transport), knowing ways of exporting the MC-results, in this case the simulated gravity and the lithology-blocks, comes in handy. There are many different ways of saving stuff (e.g. pickle the simulation results), but here we present simple exports as `.csv` and `.npy` files.

.. GENERATED FROM PYTHON SOURCE LINES 192-203

.. code-block:: default



    gravdf = pd.DataFrame.from_dict(grav)

    # add station coordinates to the dataframe
    gravdf["X"] = station_coordinates[:,0]
    gravdf["Y"] = station_coordinates[:,1]
    gravdf["Z"] =station_coordinates[:,2]

    gravdf.head()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Real_0</th>
          <th>Real_1</th>
          <th>Real_2</th>
          <th>Real_3</th>
          <th>Real_4</th>
          <th>Real_5</th>
          <th>Real_6</th>
          <th>Real_7</th>
          <th>Real_8</th>
          <th>Real_9</th>
          <th>X</th>
          <th>Y</th>
          <th>Z</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>-533.673251</td>
          <td>-532.673751</td>
          <td>-534.038468</td>
          <td>-533.143688</td>
          <td>-534.590898</td>
          <td>-533.483330</td>
          <td>-533.631910</td>
          <td>-534.186293</td>
          <td>-533.873075</td>
          <td>-534.097548</td>
          <td>21777.78</td>
          <td>7142.86</td>
          <td>405.17</td>
        </tr>
        <tr>
          <th>1</th>
          <td>-533.469712</td>
          <td>-531.718800</td>
          <td>-534.009147</td>
          <td>-532.039845</td>
          <td>-534.347122</td>
          <td>-533.275549</td>
          <td>-533.126723</td>
          <td>-534.074861</td>
          <td>-533.600943</td>
          <td>-534.057343</td>
          <td>22343.43</td>
          <td>9142.86</td>
          <td>455.98</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-529.378333</td>
          <td>-529.036879</td>
          <td>-528.462048</td>
          <td>-529.954048</td>
          <td>-525.608291</td>
          <td>-527.275759</td>
          <td>-528.747860</td>
          <td>-527.931698</td>
          <td>-527.359286</td>
          <td>-526.207417</td>
          <td>16686.87</td>
          <td>4285.71</td>
          <td>389.34</td>
        </tr>
        <tr>
          <th>3</th>
          <td>-533.691969</td>
          <td>-531.952383</td>
          <td>-534.055335</td>
          <td>-532.632034</td>
          <td>-533.845509</td>
          <td>-532.641484</td>
          <td>-533.093253</td>
          <td>-534.387168</td>
          <td>-533.574950</td>
          <td>-533.370863</td>
          <td>21494.95</td>
          <td>8000.00</td>
          <td>424.80</td>
        </tr>
        <tr>
          <th>4</th>
          <td>-533.735260</td>
          <td>-532.092981</td>
          <td>-534.197462</td>
          <td>-532.649674</td>
          <td>-533.692935</td>
          <td>-532.890857</td>
          <td>-533.219227</td>
          <td>-534.262520</td>
          <td>-533.789781</td>
          <td>-533.422763</td>
          <td>21494.95</td>
          <td>11428.57</td>
          <td>469.85</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 204-205

This can be saved as usual with `df.to_csv('pathname')` using Pandas. For the lithological block model, one good option is to save it as a numpy array, using `numpy.save()`.

.. GENERATED FROM PYTHON SOURCE LINES 205-208

.. code-block:: default


    np.save('../../data/outputs/MCexample_10realizations.npy', lith_blocks)








.. GENERATED FROM PYTHON SOURCE LINES 209-212

Quick model analysis
--------------------
Let's have a quick first look at the resulting gravity and lithological block models. From the gravity dictionary, we can quickly generate a dataframe, convenient for further model analysis.

.. GENERATED FROM PYTHON SOURCE LINES 212-217

.. code-block:: default



    prob_block = gp.bayesian.fields.probability(lith_blocks)
    ie_block = gp.bayesian.fields.information_entropy(prob_block)








.. GENERATED FROM PYTHON SOURCE LINES 218-219

The following plot shows the probability of unit 5 in the probability block. With faults not being excluded, and counting of units starting with 0, we can see that the index 5 relates to the `Lower-filling` surface. The plot shows where to expect the unit. Everywhere, this unit is present throughout the simulations, the probability plot shows a bright yellow (probability = 1). Where it is always absent, we see the dark violet (probability = 0). The blueish-greenish areas are in between, meaning that in some realizations, the `Lower-filling` unit is present there, in other realization it is not.

.. GENERATED FROM PYTHON SOURCE LINES 219-228

.. code-block:: default


    layer = 5
    gp.plot_2d(geo_model,
                show_lith=False, show_boundaries=False, show_data=False,
                regular_grid=prob_block[layer],
                kwargs_regular_grid={'cmap': 'viridis',
                                     'norm': None}
                );




.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_003.png
    :alt: Cell Number: mid Direction: y
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none


    <gempy.plot.visualization_2d.Plot2D object at 0x0000017586724700>



.. GENERATED FROM PYTHON SOURCE LINES 229-232

In the for-loop above, we not only varied the bottom boundary of the `Lower-filling` unit, but also `Upper-filling` and `Unconformity`. Using the measure of information entropy, we can visualize the parts of the model, where the most change is happening, i.e. where entropy is largest. Black areas in the following plot have zero information entropy, as there is only one "microstate" for the system, i.e. the model ensemble.  

This means, we'd always encounter the same unit at the same place in every ensemble member. The colored areas, however, are areas where we'd encounter different geological units between ensemble members.

.. GENERATED FROM PYTHON SOURCE LINES 232-240

.. code-block:: default


    gp.plot_2d(geo_model,
                show_lith=False, show_boundaries=False, show_data=False,
                regular_grid=ie_block,
                kwargs_regular_grid={'cmap': 'magma',
                                     'norm': None}
                );




.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_004.png
    :alt: Cell Number: mid Direction: y
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none


    <gempy.plot.visualization_2d.Plot2D object at 0x0000017581AF3F10>



.. GENERATED FROM PYTHON SOURCE LINES 241-242

Finally, let's have a look at the gravity. We'll simply have a look at mean and standard deviation of the simulated gravity of the ensemble:

.. GENERATED FROM PYTHON SOURCE LINES 242-269

.. code-block:: default


    # make subplots with mean and std
    gravdf_plt = pd.DataFrame.from_dict(grav)
    gravdf_plt.to_csv('../../data/outputs/MCexample_10grav.csv', index=False)

    fig, axs = plt.subplots(1,2, figsize=[15,5], sharey=True)
    m_grav = np.mean(gravdf_plt, axis=1)
    st_grav = np.std(gravdf_plt, axis=1)

    m = axs[0].scatter(grav_stations['X'], grav_stations['Y'], c=m_grav, 
               marker='s', s=90, cmap='magma', zorder=2)
    axs[0].contourf(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],20, cmap='gist_earth', zorder=0)
    axs[0].contour(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],10, colors='gray', zorder=1)
    s = axs[1].scatter(grav_stations['X'], grav_stations['Y'], c=st_grav,
                  marker='s', s=90, cmap='magma', zorder=2)
    axs[1].contourf(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],20, cmap='gist_earth', zorder=0)
    axs[1].contour(dtm[:,:,0], dtm[:,:,1], dtm[:,:,2],10, colors='gray', zorder=1)
    fig.colorbar(m, ax=axs[0], label='gravity')
    fig.colorbar(s, ax=axs[1], label='std of gravity')
    axs[0].set_title('Ensemble mean')
    axs[1].set_title('Ensemble standard deviation')
    axs[0].set_ylabel('Y [m]')
    axs[0].set_xlabel('X [m]')
    axs[1].set_xlabel('X [m]')

    fig.tight_layout()




.. image:: /Tutorials/images/sphx_glr_02_POC_create-MC-ensemble_005.png
    :alt: Ensemble mean, Ensemble standard deviation
    :class: sphx-glr-single-img






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

   **Total running time of the script:** ( 4 minutes  54.026 seconds)


.. _sphx_glr_download_Tutorials_02_POC_create-MC-ensemble.py:


.. only :: html

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



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

     :download:`Download Python source code: 02_POC_create-MC-ensemble.py <02_POC_create-MC-ensemble.py>`



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

     :download:`Download Jupyter notebook: 02_POC_create-MC-ensemble.ipynb <02_POC_create-MC-ensemble.ipynb>`


.. only:: html

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

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