Merge pull request #186 from Uncertainty-Group-Braunschweig/Development

New Sensitivity Indices

docs/code/sensitivity/chatterjee/README.rst

@@ -0,0 +1,15 @@
+Chatterjee indices
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+These examples serve as a guide for using the Chatterjee sensitivity module. They have been taken from various papers to enable validation of the implementation and have been referenced accordingly.
+
+1. **Ishigami function**
+
+    In addition to the Pick and Freeze scheme, the Sobol indices can be estimated using the rank statistics approach :cite:`gamboa2020global`. We demonstrate this estimation of the Sobol indices using the Ishigami function.
+
+2. **Exponential function**
+
+    For the Exponential model, analytical Cramér-von Mises indices are available :cite:`CVM` and since they are equivalent to the Chatterjee indices in the sample limit, they are shown here.
+
+3. **Sobol function**
+
+    This example was considered in :cite:`gamboa2020global` (page 18) to compare the Pick and Freeze scheme with the rank statistics approach for estimating the Sobol indices.

docs/code/sensitivity/chatterjee/local_exponential.py

@@ -0,0 +1,20 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+
+
+def evaluate(X: np.array) -> np.array:
+    r"""A non-linear function that is used to demonstrate sensitivity index.
+
+    .. math::
+        f(x) = \exp(x_1 + 2*x_2)
+    """
+
+    Y = np.exp(X[:, 0] + 2 * X[:, 1])
+
+    return Y

docs/code/sensitivity/chatterjee/local_ishigami.py

@@ -0,0 +1,23 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+
+
+def evaluate(X, params=[7, 0.1]):
+    """Non-monotonic Ishigami-Homma three parameter test function"""
+
+    a = params[0]
+    b = params[1]
+
+    Y = (
+        np.sin(X[:, 0])
+        + a * np.power(np.sin(X[:, 1]), 2)
+        + b * np.power(X[:, 2], 4) * np.sin(X[:, 0])
+    )
+
+    return Y

docs/code/sensitivity/chatterjee/local_sobol_func.py

@@ -0,0 +1,42 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+import copy
+
+
+def evaluate(X, a_values):
+
+    dims = len(a_values)
+    g = 1
+
+    for i in range(dims):
+        g_i = (np.abs(4 * X[:, i] - 2) + a_values[i]) / (1 + a_values[i])
+        g *= g_i
+
+    return g
+
+
+def sensitivities(a_values):
+
+    dims = len(a_values)
+
+    Total_order = np.zeros((dims, 1))
+
+    V_i = 1 / (3 * (1 + a_values) ** 2)
+
+    total_variance = np.prod(1 + V_i) - 1
+
+    First_order = V_i / total_variance
+
+    for i in range(dims):
+
+        rem_First_order = copy.deepcopy(V_i)
+        rem_First_order[i] = 0
+        Total_order[i] = V_i[i] * np.prod(rem_First_order + 1) / total_variance
+
+    return First_order.reshape(-1, 1), Total_order

docs/code/sensitivity/chatterjee/plot_chatterjee_exponential.py

@@ -0,0 +1,76 @@
+"""
+
+Exponential function
+==============================================
+
+The exponential function was used in [1]_ to demonstrate the 
+Cramér-von Mises indices. Chattererjee indices approach the Cramér-von Mises 
+indices in the sample limit and will be demonstrated via this example.
+
+.. math::
+    f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
+
+.. [1] Gamboa, F., Klein, T., & Lagnoux, A. (2018). Sensitivity Analysis Based on \
+Cramér-von Mises Distance. SIAM/ASA Journal on Uncertainty Quantification, 6(2), \
+522-548. doi:10.1137/15M1025621. (`Link <https://doi.org/10.1137/15M1025621>`_)
+
+"""
+
+# %%
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Normal
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.Chatterjee import Chatterjee
+from UQpy.sensitivity.PostProcess import *
+
+np.random.seed(123)
+
+# %% [markdown]
+# **Define the model and input distributions**
+
+# Create Model object
+model = PythonModel(
+    model_script="local_exponential.py",
+    model_object_name="evaluate",
+    var_names=[
+        "X_1",
+        "X_2",
+    ],
+    delete_files=True,
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Normal(0, 1)] * 2)
+
+# %% [markdown]
+# **Compute Chatterjee indices**
+
+# %% [markdown]
+SA = Chatterjee(runmodel_obj, dist_object)
+
+# Compute Chatterjee indices using the pick and freeze algorithm
+computed_indices = SA.run(n_samples=1_000_000)
+
+# %% [markdown]
+# **Chattererjee indices**
+#
+# Chattererjee indices approach the Cramér-von Mises indices in the sample limit.
+#
+# Expected value of the sensitivity indices:
+#
+# :math:`S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145`
+#
+# :math:`S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693`
+
+# %%
+computed_indices["chatterjee_i"]
+
+# **Plot the Chatterjee indices**
+fig1, ax1 = plot_sensitivity_index(
+    computed_indices["chatterjee_i"][:, 0],
+    plot_title="Chatterjee indices",
+    color="C2",
+)

docs/code/sensitivity/chatterjee/plot_chatterjee_ishigami.py

@@ -0,0 +1,98 @@
+r"""
+
+Ishigami function
+==============================================
+
+The ishigami function is a non-linear, non-monotonic function that is commonly used to 
+benchmark uncertainty and senstivity analysis methods.
+
+.. math::
+    f(x_1, x_2, x_3) = sin(x_1) + a \cdot sin^2(x_2) + b \cdot x_3^4 sin(x_1)
+
+.. math::
+    x_1, x_2, x_3 \sim \mathcal{U}(-\pi, \pi), \quad a, b\in \mathbb{R}
+
+"""
+
+# %%
+import numpy as np
+
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Uniform
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.Chatterjee import Chatterjee
+from UQpy.sensitivity.PostProcess import *
+
+np.random.seed(123)
+
+# %% [markdown]
+# **Define the model and input distributions**
+
+# Create Model object
+model = PythonModel(
+    model_script="local_ishigami.py",
+    model_object_name="evaluate",
+    var_names=[r"$X_1$", "$X_2$", "$X_3$"],
+    delete_files=True,
+    params=[7, 0.1],
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Uniform(-np.pi, 2 * np.pi)] * 3)
+
+# %% [markdown]
+# **Compute Chatterjee indices**
+
+# %% [markdown]
+SA = Chatterjee(runmodel_obj, dist_object)
+
+computed_indices = SA.run(
+    n_samples=100_000,
+    estimate_sobol_indices=True,
+    num_bootstrap_samples=100,
+    confidence_level=0.95,
+)
+
+# %% [markdown]
+# **Chattererjee indices**
+
+# %%
+computed_indices["chatterjee_i"]
+
+# %% [markdown]
+# **Confidence intervals for the Chatterjee indices**
+
+# %%
+computed_indices["confidence_interval_chatterjee_i"]
+
+# **Plot the Chatterjee indices**
+fig1, ax1 = plot_sensitivity_index(
+    computed_indices["chatterjee_i"][:, 0],
+    computed_indices["confidence_interval_chatterjee_i"],
+    plot_title="Chatterjee indices",
+    color="C2",
+)
+
+# %% [markdown]
+# **Estimated Sobol indices**
+#
+# Expected first order Sobol indices:
+#
+# :math:`S_1`: 0.3139
+#
+# :math:`S_2`: 0.4424
+#
+# :math:`S_3`: 0.0
+
+# %%
+computed_indices["sobol_i"]
+
+# **Plot the first order Sobol indices**
+fig2, ax2 = plot_sensitivity_index(
+    computed_indices["sobol_i"][:, 0],
+    plot_title="First order Sobol indices",
+    color="C0",
+)