Skip to content

Commit

Permalink
Added examples Chatterjee sensitivity
Browse files Browse the repository at this point in the history
  • Loading branch information
Prateek Bhustali committed May 8, 2022
1 parent 0b5a666 commit e4407c6
Show file tree
Hide file tree
Showing 6 changed files with 267 additions and 0 deletions.
20 changes: 20 additions & 0 deletions docs/code/sensitivity/chatterjee/local_exponential.py
Original file line number Diff line number Diff line change
@@ -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
23 changes: 23 additions & 0 deletions docs/code/sensitivity/chatterjee/local_ishigami.py
Original file line number Diff line number Diff line change
@@ -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
42 changes: 42 additions & 0 deletions docs/code/sensitivity/chatterjee/local_sobol_func.py
Original file line number Diff line number Diff line change
@@ -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
54 changes: 54 additions & 0 deletions docs/code/sensitivity/chatterjee/plot_chatterjee_exponential.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
"""
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.chatterjee import Chatterjee

# %%
# 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

# %%
SA = Chatterjee(runmodel_obj, dist_object)

# Compute Sobol indices using the pick and freeze algorithm
computed_indices = SA.run(n_samples=1_000_000)

# %% [markdown]
# Cramer-von Mises sensitivity analysis
#
# Expected value of the sensitivity indices:
#
# $S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145$
#
# $S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693$

# %%
computed_indices["chatterjee_i"]
58 changes: 58 additions & 0 deletions docs/code/sensitivity/chatterjee/plot_chatterjee_ishigami.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
r"""
Ishigami function
==============================================
.. 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

# %%
# 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

# %%
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,
)

# %%
computed_indices["chatterjee_i"]

# %%
computed_indices["CI_chatterjee_i"]

# %%
computed_indices["sobol_i"]
70 changes: 70 additions & 0 deletions docs/code/sensitivity/chatterjee/plot_chatterjee_sobol_func.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
r"""
Sobol function
==============================================
.. math::
g(x_1, x_2, \ldots, x_D) := \prod_{i=1}^{D} \frac{|4x_i - 2| + a_i}{1 + a_i},
where,
.. math::
x_i \sim \mathcal{U}(0, 1), \quad a_i \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

# %%
# Create Model object
num_vars = 6
a_vals = np.array([0.0, 0.5, 3.0, 9.0, 99.0, 99.0])

model = PythonModel(
model_script="local_sobol_func.py",
model_object_name="evaluate",
var_names=["X_" + str(i) for i in range(num_vars)],
delete_files=True,
a_values=a_vals,
)

runmodel_obj = RunModel(model=model)

# Define distribution object
dist_object = JointIndependent([Uniform(0, 1)] * num_vars)

# %% [markdown]
# Compute Chatterjee indices

# %%
SA = Chatterjee(runmodel_obj, dist_object)

# Compute Sobol indices using the pick and freeze algorithm
computed_indices = SA.run(n_samples=500_000, estimate_sobol_indices=True)

# %%
computed_indices["chatterjee_i"]

# %% [markdown]
# $S_1$ = 5.86781190e-01
#
# $S_2$ = 2.60791640e-01
#
# $S_3$ = 3.66738244e-02
#
# $S_4$ = 5.86781190e-03
#
# $S_5$ = 5.86781190e-05
#
# $S_6$ = 5.86781190e-05

# %%
computed_indices["sobol_i"]

0 comments on commit e4407c6

Please sign in to comment.