Creating a new stochastic example

import xarray as xr

input_data = xr.Dataset({

    "thickness_mean": ((), 300),

    "thickness_sd": ((), 50),

    "ntg": ((), 0.5),

    "porosity": ((), 0.5),

    "depth": ((), 5000),

    "ln_permeability_mean": ((), 5),

    "ln_permeability_sd": ((), 0.5),

    "thickness_mean": 300,

    "thickness_sd": 50,

    "ntg": 0.5,

    "porosity": 0.5,

    "depth": 5000,

    "ln_permeability_mean": 5,

    "ln_permeability_sd": 0.5,

})

utc_properties = instantiate_utc_properties_from_xml("path/to/valid/xml/file")

##  Portfolio run and analysis of doublet simulations

This is simple example of how to use the `calculate_doublet_performance_stochastic` function from 

the `pythermogis` package to run the simulation for a set of maps and

p-values. 

the `pythermogis` package to run the simulation for a set of maps and p-values. 

This example demonstrates the following steps:

    plt.tight_layout()  # ensure there is enough spacing

    plt.savefig(output_data_path / "npv.png")

def test_example8():

    from pythermogis import calculate_doublet_performance

    import xarray as xr

    from matplotlib import pyplot as plt

    # generate simulation samples across desired reservoir properties

    Nsamples = 1000

    thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)

    porosity_samples = np.random.uniform(low=0.1, high=0.8, size=Nsamples)

    ntg_samples = np.random.uniform(low=0.1, high=0.5, size=Nsamples)

    depth_samples = np.random.normal(loc=2000, scale=100, size=Nsamples)

    permeability_samples = np.random.normal(loc=400, scale=100, size=Nsamples)

    reservoir_properties = xr.Dataset(

        {

            "thickness": (["sample"], thickness_samples),

            "porosity": (["sample"], porosity_samples),

            "ntg": (["sample"], ntg_samples),

            "depth": (["sample"], depth_samples),

            "permeability": (["sample"], permeability_samples),

        },

        coords={"sample": np.arange(Nsamples)}

    )

    results = calculate_doublet_performance(reservoir_properties, print_execution_duration=True)

    # plot results against permeability

    fig, axes = plt.subplots(2,2, figsize=(10,10))

    plt.suptitle("Simulation Samples")

    results.plot.scatter(x="permeability", y="transmissivity", ax=axes[0,0])

    results.plot.scatter(x="permeability", y="utc", ax=axes[0,1])

    results.plot.scatter(x="permeability", y="power", ax=axes[1,0])

    results.plot.scatter(x="permeability", y="npv", ax=axes[1,1])

    plt.show()

    # Calculate Stochastics Across Samples

    def plot_survival_function(samples, xlabel='Parameter value (x)', ax=None):

        sorted_samples = np.sort(samples)

        cdf = np.arange(1, len(samples) + 1) / len(samples)

        # Compute P(θ > x)

        survival = 1 - cdf

        # Plot survival function (expectation curve as P(>x))

        if ax is None:

            fig, ax = plt.subplots(figsize=(8, 5))

        ax.plot(sorted_samples, survival, label=r'$P(\theta > x)$', color='blue')

        ax.set_xlabel(xlabel)

        ax.set_ylabel('Cumulative probability $P(\\theta > x)$')

        ax.grid(True)

        ax.legend()

    fig, axes = plt.subplots(2,2, figsize=(10,10))

    plt.suptitle("Stochastic Measures")

    plot_survival_function(results.permeability, xlabel="permeability", ax=axes[0,0])

    plot_survival_function(results.power, xlabel="power", ax=axes[0, 1])

    plot_survival_function(results.utc, xlabel="utc", ax=axes[1, 0])

    plot_survival_function(results.npv, xlabel="net-present-value", ax=axes[1, 1])

    plt.show()