Added examples for Sobol indices

docs/code/sensitivity/sobol/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/sobol/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/sobol/local_mechanical_oscillator_ODE.py

@@ -0,0 +1,60 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+from scipy.integrate import solve_ivp
+
+
+def mech_oscillator(input_parameters):
+    """
+    We have the second order differential equation:
+
+    .. math::
+
+        m \ddot{x} + c \dot{x} + k x = 0
+
+    with initial conditions: :math: `x(0) = \ell`, :math: `\dot{x}(0) = 0`.
+
+    where, for example  :math: `m \sim \mathcal{U}(10, 12)`,
+                        :math: `c \sim \mathcal{U}(0.4, 0.8)`
+                        :math: `k \sim \mathcal{U}(70, 90)`
+                        :math: `\ell \sim \mathcal{U}(-1, -0.25)`.
+
+
+    References
+    ----------
+
+    .. [1] Gamboa F, Janon A, Klein T, Lagnoux A, others .
+        Sensitivity analysis for multidimensional and functional outputs.
+        Electronic journal of statistics 2014; 8(1): 575-603.
+
+    """
+
+    # unpack the input parameters
+    m, c, k, l = input_parameters[0]
+
+    # intial conditions
+    x_0 = l
+    v_0 = 0
+
+    # time points
+    t_0 = 0
+    t_f = 40
+    dt = 0.05
+    n_t = int((t_f - t_0) / dt)
+    T = np.linspace(t_0, t_f, n_t)
+
+    def ODE(t, y):
+        """
+        The ODE system.
+        """
+        return np.array([y[1], -(k / m) * y[0] - (c / m) * y[1]])
+
+    # solve the ODE
+    sol = solve_ivp(ODE, [t_0, t_f], [x_0, v_0], method="RK45", t_eval=T)
+
+    return sol.y[0]

docs/code/sensitivity/sobol/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 = (3 * (1 + a_values) ** 2) ** (-1)
+
+    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/sobol/plot_mechanical_oscillator_ODE.py

@@ -0,0 +1,92 @@
+r"""
+
+Mechanical oscillator model (multioutput)
+==============================================
+
+The mechanical oscillator is governed by the following second-order ODE:
+
+.. math::
+    m \ddot{x} + c \dot{x} + k x = 0
+
+.. math::
+    x(0) = \ell, \dot{x}(0) = 0.
+
+The parameteres of the oscillator are modeled as follows:
+
+.. math::
+    m \sim \mathcal{U}(10, 12), c \sim \mathcal{U}(0.4, 0.8), k \sim \mathcal{U}(70, 90), \ell \sim \mathcal{U}(-1, -0.25).
+
+"""
+
+# %%
+import numpy as np
+import matplotlib.pyplot as plt
+
+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.sobol import Sobol
+
+# %%
+# Create Model object
+model = PythonModel(
+    model_script="local_mechanical_oscillator_ODE.py",
+    model_object_name="mech_oscillator",
+    var_names=[r"$m$", "$c$", "$k$", "$\ell$"],
+    delete_files=True,
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+M = Uniform(10, (12 - 10))
+C = Uniform(0.4, (0.8 - 0.4))
+K = Uniform(70, (90 - 70))
+L = Uniform(-1, (-0.25 - -1))
+dist_object = JointIndependent([M, C, K, L])
+
+# %%
+SA = Sobol(runmodel_obj, dist_object)
+
+computed_indices = SA.run(n_samples=500)
+
+# %%
+# Plot the Sobol indices
+t_0 = 0
+t_f = 40
+dt = 0.05
+n_t = int((t_f - t_0) / dt)
+T = np.linspace(t_0, t_f, n_t)
+
+fig, ax = plt.subplots(1, 2, figsize=(16, 8))
+
+ax[0].plot(T, computed_indices["sobol_total_i"][0, :], "r", label=r"$m$")
+ax[0].plot(T, computed_indices["sobol_total_i"][1, :], "g", label=r"$c$")
+ax[0].plot(T, computed_indices["sobol_total_i"][2, :], label=r"$k$", color="royalblue")
+ax[0].plot(
+    T, computed_indices["sobol_total_i"][3, :], label=r"$\ell$", color="aquamarine"
+)
+
+ax[0].set_title("Total order Sobol indices", fontsize=16)
+ax[0].set_xlabel("time (s)", fontsize=16)
+ax[0].set_ylabel(r"$S_{T_i}$", fontsize=16)
+ax[0].set_xbound(0, t_f)
+ax[0].set_ybound(-0.2, 1.2)
+ax[0].legend()
+
+ax[1].plot(T, computed_indices["sobol_i"][0, :], "r", label=r"$m$")
+ax[1].plot(T, computed_indices["sobol_i"][1, :], "g", label=r"$c$")
+ax[1].plot(T, computed_indices["sobol_i"][2, :], label=r"$k$", color="royalblue")
+ax[1].plot(T, computed_indices["sobol_i"][3, :], label=r"$\ell$", color="aquamarine")
+
+ax[1].set_title("First order Sobol indices", fontsize=16)
+ax[1].set_xlabel("time (s)", fontsize=16)
+ax[1].set_ylabel(r"$S_i$", fontsize=16)
+ax[1].set_xbound(0, t_f)
+ax[1].set_ybound(-0.2, 1.2)
+ax[1].legend(fontsize=12)
+
+fig.suptitle("Pointwise-in-time Sobol indices", fontsize=20)
+
+plt.show()

docs/code/sensitivity/sobol/plot_sobol_exponential.py

@@ -0,0 +1,60 @@
+"""
+
+Exponential function
+==============================================
+
+.. math::
+    f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
+
+"""
+
+# %%
+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.sobol import Sobol
+
+# %%
+# 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 Sobol indices
+
+# %%
+SA = Sobol(runmodel_obj, dist_object)
+
+# Compute Sobol indices using the pick and freeze algorithm
+computed_indices = SA.run(
+    n_samples=100_000, num_bootstrap_samples=1_000, confidence_level=0.95
+)
+
+# %% [markdown]
+# Expected first order Sobol indices (computed analytically):
+#
+# X1: 0.0118
+#
+# X2: 0.3738
+
+# %%
+computed_indices["sobol_i"]
+
+# %% [markdown]
+# Confidence intervals for first order Sobol indices
+
+# %%
+computed_indices["CI_sobol_i"]