From 29b47271e500ed76dfbab00156831536957be1d0 Mon Sep 17 00:00:00 2001 From: Prateek Bhustali Date: Mon, 9 May 2022 03:24:53 +0200 Subject: [PATCH] Add unit tests for generalised sobol sensitivity --- tests/unit_tests/sensitivity/multioutput.py | 42 +++ .../sensitivity/test_generalised_sobol.py | 315 ++++++++++++++++++ 2 files changed, 357 insertions(+) create mode 100644 tests/unit_tests/sensitivity/multioutput.py create mode 100644 tests/unit_tests/sensitivity/test_generalised_sobol.py diff --git a/tests/unit_tests/sensitivity/multioutput.py b/tests/unit_tests/sensitivity/multioutput.py new file mode 100644 index 00000000..6974d37c --- /dev/null +++ b/tests/unit_tests/sensitivity/multioutput.py @@ -0,0 +1,42 @@ +"""" +This is the toy example with multiple outputs from [1]_. + +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. + +""" + +import numpy as np + + +def evaluate(X): + + """ + + * **Input:** + + * **X** (`ndarray`): + Samples from the input distribution. + Shape: (n_samples, 2) + + * **Output:** + + * **Y** (`ndarray`): + Model evaluations. + Shape: (2, n_samples) + + """ + + n_samples = X.shape[0] + + output = np.zeros((2, n_samples)) + + output[0, :] = X[:, 0] + X[:, 1] + X[:, 0] * X[:, 1] + + output[1, :] = 2 * X[:, 0] + X[:, 1] + 3 * X[:, 0] * X[:, 1] + + return output diff --git a/tests/unit_tests/sensitivity/test_generalised_sobol.py b/tests/unit_tests/sensitivity/test_generalised_sobol.py new file mode 100644 index 00000000..c759d85b --- /dev/null +++ b/tests/unit_tests/sensitivity/test_generalised_sobol.py @@ -0,0 +1,315 @@ +"""" +This is the test module for the Generalised Sobol indices. + +Here, we will use the toy example from [1]_, which is a multi-output problem. + + +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. + +Important +---------- +The computed indices are computed using the `np.isclose` function. + +Function signature: + numpy.isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False) + + Parameters: + a, b: array_like + Input arrays to compare. + + rtol: float + The relative tolerance parameter. + + atol: float + The absolute tolerance parameter. + +Each element of the `diff` array is compared as follows: +diff = |a - b| +diff <= atol + rtol * abs(b) + +- relative tolerance: rtol * abs(b) + It is the maximum allowed difference between a and b, + relative to the absolute value of b. + For example, to set a tolerance of 1%, pass rol=0.01, + which assures that the values are within 2 decimal places of each other. +- absolute tolerance: atol + When b is close to zero, the atol value is used. + +""" + +import numpy as np +import pytest +import scipy + +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel +from UQpy.distributions import Uniform, Normal +from UQpy.distributions.collection.JointIndependent import JointIndependent +from UQpy.sensitivity.generalised_sobol import GeneralisedSobol + +# Prepare +############################################################################### + +# Prepare the input distribution +@pytest.fixture() +def normal_input_dist_object(): + """ + This function returns the input distribution for the toy model. + + X1 ~ Normal(0, 1) + X2 ~ Normal(0, 1) + + """ + return JointIndependent([Normal(0, 1)] * 2) + + +@pytest.fixture() +def uniform_input_dist_object(): + """ + This function returns the input distribution for the toy model. + + X1 ~ Uniform(0, 1) + X2 ~ Uniform(0, 1) + + """ + return JointIndependent([Uniform(0, 1)] * 2) + + +@pytest.fixture() +def toy_model_object(): + """ + This function creates the toy model. + + """ + model = PythonModel( + model_script="multioutput.py", + model_object_name="evaluate", + var_names=[ + "X_1", + "X_2", + ], + delete_files=True, + ) + + runmodel_obj = RunModel(model=model) + + return runmodel_obj + + +@pytest.fixture() +def generalised_sobol_object_normal(normal_input_dist_object, toy_model_object): + """ + This function creates the Generalised Sobol indices object + with normal input distribution. + + """ + + return GeneralisedSobol(toy_model_object, normal_input_dist_object) + + +@pytest.fixture() +def generalised_sobol_object_uniform(uniform_input_dist_object, toy_model_object): + """ + This function creates the Generalised Sobol indices object + with uniform input distribution. + + """ + + return GeneralisedSobol(toy_model_object, uniform_input_dist_object) + + +@pytest.fixture() +def analytical_toy_GSI_normal(): + """ + Analytical first order Generalised Sobol indices + for the toy example with normal input distribution. + """ + + return np.array([0.2941, 0.1176]).reshape(-1, 1) + + +@pytest.fixture() +def analytical_toy_GSI_uniform(): + """ " + Analytical first order Generalised Sobol indices + for toy example with uniform input distribution. + """ + + return np.array([0.6084, 0.3566]).reshape(-1, 1) + + +@pytest.fixture() +def pick_and_freeze_toy_GSI_normal(generalised_sobol_object_normal): + """ " + Generalised first order Sobol indices computed using the Pick and Freeze + approach for the toy example with normal input distribution. + """ + + SA = generalised_sobol_object_normal + + np.random.seed(12345) #! set seed for reproducibility + + computed_indices = SA.run(n_samples=100_000) + + return computed_indices["gen_sobol_i"] + + +@pytest.fixture() +def pick_and_freeze_toy_GSI_uniform(generalised_sobol_object_uniform): + """ " + Generalised first order Sobol indices computed using the Pick and Freeze + approach for the toy example with uniform input distribution. + """ + + SA = generalised_sobol_object_uniform + + np.random.seed(12345) #! set seed for reproducibility + + computed_indices = SA.run(n_samples=100_000) + + return computed_indices["gen_sobol_i"] + + +@pytest.fixture() +def NUM_SAMPLES(): + """This function returns the number of samples for bootstrapping""" + + num_bootstrap_samples = 500 + num_samples = 20_000 + + return num_bootstrap_samples, num_samples + + +@pytest.fixture() +def bootstrap_generalised_sobol_index_variance( + generalised_sobol_object_normal, NUM_SAMPLES +): + + SA = generalised_sobol_object_normal + + np.random.seed(12345) #! set seed for reproducibility + + num_bootstrap_samples, n_samples = NUM_SAMPLES + + confidence_level = 0.95 + delta = -scipy.stats.norm.ppf((1 - confidence_level) / 2) + + # Compute the confidence intervals + + computed_indices = SA.run( + n_samples=n_samples, + num_bootstrap_samples=num_bootstrap_samples, + confidence_level=confidence_level, + ) + + gen_sobol_i = computed_indices["gen_sobol_i"].ravel() + gen_sobol_total_i = computed_indices["gen_sobol_total_i"].ravel() + upper_bound_first_order = computed_indices["CI_gen_sobol_i"][:, 1] + upper_bound_total_order = computed_indices["CI_gen_sobol_total_i"][:, 1] + + std_bootstrap_first_order = (upper_bound_first_order - gen_sobol_i) / delta + std_bootstrap_total_order = (upper_bound_total_order - gen_sobol_total_i) / delta + + return std_bootstrap_first_order**2, std_bootstrap_total_order**2 + + +@pytest.fixture() +def model_eval_generalised_sobol_index_variance(): + + """ + For computational efficiency, the variance of the generalised Sobol indices + is precomputed using model evaluations with + NUM_SAMPLES (num_repetitions=500, num_samples=20_000) + + Copy-paste the following code to generate the variance + of the Sobol indices: + + runmodel_obj = RunModel(model_script='multioutput.py', + model_object_name='multioutput_toy', + vec=True, delete_files=True) + + dist_object_1 = JointInd([Normal(0, 1)]*2) + + SA = GeneralisedSobol(runmodel_obj, dist_object_1) + + np.random.seed(12345) # for reproducibility + + num_repetitions, n_samples = 500, 20_000 + + num_vars = 2 + + bootstrap_first_order = np.zeros((num_vars, num_bootstrap_samples)) + bootstrap_total_order = np.zeros((num_vars, num_bootstrap_samples)) + + for b in range(num_repetitions): + + computed_indices = SA.run(n_samples=n_samples) + + bootstrap_first_order[:, b] = computed_indices["gen_sobol_i"].ravel() + bootstrap_total_order[:, b] = computed_indices["gen_sobol_total_i"].ravel() + + var_bootstrap_gen_S = np.var(bootstrap_first_order, axis=1, ddof=1) + var_bootstrap_gen_S_T = np.var(bootstrap_total_order, axis=1, ddof=1) + + print(var_bootstrap_gen_S) + print(var_bootstrap_gen_S_T) + + """ + + variance_first_order = np.array([0.00011284, 0.00012608]) + + variance_total_order = np.array([0.00012448, 0.00011208]) + + return variance_first_order, variance_total_order + + +# Unit tests +############################################################################### + + +def test_pick_and_freeze_estimator( + pick_and_freeze_toy_GSI_normal, + analytical_toy_GSI_normal, + pick_and_freeze_toy_GSI_uniform, + analytical_toy_GSI_uniform, +): + """ + Test the pick and freeze estimator. + + """ + + # Prepare + N_true = analytical_toy_GSI_normal + N_estimate = pick_and_freeze_toy_GSI_normal + + U_true = analytical_toy_GSI_uniform + U_estimate = pick_and_freeze_toy_GSI_uniform + + # Act + # Idea: Measure accuracy upto 2 decimal places -> rtol=0, atol=1e-2 + assert np.isclose(N_estimate, N_true, rtol=0, atol=1e-2).all() + assert np.isclose(U_estimate, U_true, rtol=0, atol=1e-2).all() + + +def test_bootstrap_variance_computation( + model_eval_generalised_sobol_index_variance, + bootstrap_generalised_sobol_index_variance, +): + + """Test the bootstrap variance computation.""" + + # Prepare + var_first, var_total = model_eval_generalised_sobol_index_variance + boot_var_first, boot_var_total = bootstrap_generalised_sobol_index_variance + + # Act + assert var_first.shape == boot_var_first.shape + + # Idea: Ensure bootstrap variance and MC variance are of same order -> rtol=0, atol=1e-4 + assert np.isclose(boot_var_first, var_first, rtol=0, atol=1e-4).all() + assert np.isclose(boot_var_total, var_total, rtol=0, atol=1e-4).all()