Morris analysis

from __future__ import annotations

import pprint

from gemseo import configure_logger
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

configure_logger()

<RootLogger root (INFO)>

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()

 WARNING - 08:35:12: No coupling in MDA, switching chain_linearize to True.
    INFO - 08:35:12:
    INFO - 08:35:12: *** Start MorrisAnalysisSamplingPhase execution ***
    INFO - 08:35:12: MorrisAnalysisSamplingPhase
    INFO - 08:35:12:    Disciplines: IshigamiDiscipline
    INFO - 08:35:12:    MDO formulation: MDF
    INFO - 08:35:12: Running the algorithm MorrisDOE:
    INFO - 08:35:12:      5%|▌         | 1/20 [00:00<00:00, 367.89 it/sec]
    INFO - 08:35:12:     10%|█         | 2/20 [00:00<00:00, 585.47 it/sec]
    INFO - 08:35:12:     15%|█▌        | 3/20 [00:00<00:00, 745.70 it/sec]
    INFO - 08:35:12:     20%|██        | 4/20 [00:00<00:00, 860.11 it/sec]
    INFO - 08:35:12:     25%|██▌       | 5/20 [00:00<00:00, 957.08 it/sec]
    INFO - 08:35:12:     30%|███       | 6/20 [00:00<00:00, 1036.53 it/sec]
    INFO - 08:35:12:     35%|███▌      | 7/20 [00:00<00:00, 1098.03 it/sec]
    INFO - 08:35:12:     40%|████      | 8/20 [00:00<00:00, 1140.69 it/sec]
    INFO - 08:35:12:     45%|████▌     | 9/20 [00:00<00:00, 1184.91 it/sec]
    INFO - 08:35:12:     50%|█████     | 10/20 [00:00<00:00, 1224.58 it/sec]
    INFO - 08:35:12:     55%|█████▌    | 11/20 [00:00<00:00, 1254.04 it/sec]
    INFO - 08:35:12:     60%|██████    | 12/20 [00:00<00:00, 1279.69 it/sec]
    INFO - 08:35:12:     65%|██████▌   | 13/20 [00:00<00:00, 1301.83 it/sec]
    INFO - 08:35:12:     70%|███████   | 14/20 [00:00<00:00, 1322.05 it/sec]
    INFO - 08:35:12:     75%|███████▌  | 15/20 [00:00<00:00, 1331.10 it/sec]
    INFO - 08:35:12:     80%|████████  | 16/20 [00:00<00:00, 1345.06 it/sec]
    INFO - 08:35:12:     85%|████████▌ | 17/20 [00:00<00:00, 1362.85 it/sec]
    INFO - 08:35:12:     90%|█████████ | 18/20 [00:00<00:00, 1380.21 it/sec]
    INFO - 08:35:12:     95%|█████████▌| 19/20 [00:00<00:00, 1394.53 it/sec]
    INFO - 08:35:12:    100%|██████████| 20/20 [00:00<00:00, 1405.64 it/sec]
    INFO - 08:35:12: *** End MorrisAnalysisSamplingPhase execution (time: 0:00:00.024414) ***

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