Morris analysis

from __future__ import annotations

import pprint

from gemseo.problems.uncertainty.ishigami.ishigami_discipline import IshigamiDiscipline
from gemseo.problems.uncertainty.ishigami.ishigami_space import IshigamiSpace
from gemseo.uncertainty.sensitivity.morris_analysis import MorrisAnalysis

In this example, we consider the Ishigami function [IH90]

\[f(x_1,x_2,x_3)=\sin(x_1)+7\sin(x_2)^2+0.1x_3^4\sin(x_1)\]

implemented as an Discipline by the IshigamiDiscipline. It is commonly used with the independent random variables \(X_1\), \(X_2\) and \(X_3\) uniformly distributed between \(-\pi\) and \(\pi\) and defined in the IshigamiSpace.

discipline = IshigamiDiscipline()
uncertain_space = IshigamiSpace()

Then, we run sensitivity analysis of type MorrisAnalysis:

sensitivity_analysis = MorrisAnalysis()
sensitivity_analysis.compute_samples([discipline], uncertain_space, n_samples=0)
sensitivity_analysis.compute_indices()

    INFO - 16:22:32: *** Start MorrisAnalysisSamplingPhase execution ***
    INFO - 16:22:32: MorrisAnalysisSamplingPhase
    INFO - 16:22:32:    Disciplines: IshigamiDiscipline
    INFO - 16:22:32:    MDO formulation: MDF
    INFO - 16:22:32: Running the algorithm MorrisDOE:
    INFO - 16:22:32:      5%|▌         | 1/20 [00:00<00:00, 489.02 it/sec]
    INFO - 16:22:32:     10%|█         | 2/20 [00:00<00:00, 822.82 it/sec]
    INFO - 16:22:32:     15%|█▌        | 3/20 [00:00<00:00, 1097.79 it/sec]
    INFO - 16:22:32:     20%|██        | 4/20 [00:00<00:00, 1320.21 it/sec]
    INFO - 16:22:32:     25%|██▌       | 5/20 [00:00<00:00, 1511.35 it/sec]
    INFO - 16:22:32:     30%|███       | 6/20 [00:00<00:00, 1685.25 it/sec]
    INFO - 16:22:32:     35%|███▌      | 7/20 [00:00<00:00, 1829.18 it/sec]
    INFO - 16:22:32:     40%|████      | 8/20 [00:00<00:00, 1951.52 it/sec]
    INFO - 16:22:32:     45%|████▌     | 9/20 [00:00<00:00, 2064.35 it/sec]
    INFO - 16:22:32:     50%|█████     | 10/20 [00:00<00:00, 2170.63 it/sec]
    INFO - 16:22:32:     55%|█████▌    | 11/20 [00:00<00:00, 2217.82 it/sec]
    INFO - 16:22:32:     60%|██████    | 12/20 [00:00<00:00, 2282.51 it/sec]
    INFO - 16:22:32:     65%|██████▌   | 13/20 [00:00<00:00, 2359.31 it/sec]
    INFO - 16:22:32:     70%|███████   | 14/20 [00:00<00:00, 2435.92 it/sec]
    INFO - 16:22:32:     75%|███████▌  | 15/20 [00:00<00:00, 2499.29 it/sec]
    INFO - 16:22:32:     80%|████████  | 16/20 [00:00<00:00, 2545.18 it/sec]
    INFO - 16:22:32:     85%|████████▌ | 17/20 [00:00<00:00, 2607.25 it/sec]
    INFO - 16:22:32:     90%|█████████ | 18/20 [00:00<00:00, 2665.78 it/sec]
    INFO - 16:22:32:     95%|█████████▌| 19/20 [00:00<00:00, 2711.07 it/sec]
    INFO - 16:22:32:    100%|██████████| 20/20 [00:00<00:00, 2706.44 it/sec]
    INFO - 16:22:32: *** End MorrisAnalysisSamplingPhase execution ***

MorrisAnalysis.SensitivityIndices(mu={'y': [{'x1': array([-0.60047199]), 'x2': array([0.51230435]), 'x3': array([-0.89800793])}]}, mu_star={'y': [{'x1': array([0.69887482]), 'x2': array([0.65136343]), 'x3': array([0.89805157])}]}, sigma={'y': [{'x1': array([0.96395158]), 'x2': array([0.6549141]), 'x3': array([0.79878356])}]}, relative_sigma={'y': [{'x1': array([1.37929075]), 'x2': array([1.00545113]), 'x3': array([0.88946291])}]}, min={'y': [{'x1': array([0.0338188]), 'x2': array([0.11821721]), 'x3': array([8.72820113e-05])}]}, max={'y': [{'x1': array([2.2360336]), 'x2': array([1.25769934]), 'x3': array([2.12052546])}]})

The resulting indices are the empirical means and the standard deviations of the absolute output variations due to input changes.

sensitivity_analysis.indices

MorrisAnalysis.SensitivityIndices(mu={'y': [{'x1': array([-0.60047199]), 'x2': array([0.51230435]), 'x3': array([-0.89800793])}]}, mu_star={'y': [{'x1': array([0.69887482]), 'x2': array([0.65136343]), 'x3': array([0.89805157])}]}, sigma={'y': [{'x1': array([0.96395158]), 'x2': array([0.6549141]), 'x3': array([0.79878356])}]}, relative_sigma={'y': [{'x1': array([1.37929075]), 'x2': array([1.00545113]), 'x3': array([0.88946291])}]}, min={'y': [{'x1': array([0.0338188]), 'x2': array([0.11821721]), 'x3': array([8.72820113e-05])}]}, max={'y': [{'x1': array([2.2360336]), 'x2': array([1.25769934]), 'x3': array([2.12052546])}]})

The main indices corresponds to these empirical means (this main method can be changed with MorrisAnalysis.main_method):

pprint.pprint(sensitivity_analysis.main_indices)

{'y': [{'x1': array([0.69887482]),
        'x2': array([0.65136343]),
        'x3': array([0.89805157])}]}

and can be interpreted with respect to the empirical bounds of the outputs:

pprint.pprint(sensitivity_analysis.outputs_bounds)

{'y': (array([-1.42959705]), array([14.89344259]))}

We can also get the input parameters sorted by decreasing order of influence:

sensitivity_analysis.sort_input_variables("y")

['x3', 'x1', 'x2']

We can use the method MorrisAnalysis.plot() to visualize the different series of indices:

sensitivity_analysis.plot("y", save=False, show=True, lower_mu=0, lower_sigma=0)

<Figure size 640x480 with 1 Axes>

Lastly, the sensitivity indices can be exported to a Dataset:

sensitivity_analysis.to_dataset()

GROUP	mu	mu_star	sigma	relative_sigma	min	max
VARIABLE	y	y	y	y	y	y
COMPONENT	0	0	0	0	0	0
x1	-0.600472	0.698875	0.963952	1.379291	0.033819	2.236034
x2	0.512304	0.651363	0.654914	1.005451	0.118217	1.257699
x3	-0.898008	0.898052	0.798784	0.889463	0.000087	2.120525