diff --git a/src/gemseo/uncertainty/sensitivity/fast_analysis.py b/src/gemseo/uncertainty/sensitivity/fast_analysis.py new file mode 100644 index 0000000000000000000000000000000000000000..546152660c9a17c0128a921c2b43a1a9228b7ba1 --- /dev/null +++ b/src/gemseo/uncertainty/sensitivity/fast_analysis.py @@ -0,0 +1,431 @@ +# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License version 3 as published by the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# Contributors: +# INITIAL AUTHORS - initial API and implementation and/or initial +# documentation +# :author: Gokulnath Ramachandran +# OTHER AUTHORS - MACROSCOPIC CHANGES +r"""Fourier Amplitude Sensitivity Testing (FAST). + +Notes: + ----- + FAST is a sensitivity analysis method which is based upon the ANOVA + decomposition of the variance of the model response :math:`y = f(\vect{X})`, + the latter being represented by its Fourier expansion. + :math:`\vect{X}=\{X^1,\dots,X^{n_X}\}` is an input random vector of :math:`n_X` + independent components. + + OpenTURNS implements the extended FAST method consisting in computing + alternately the first order and the total-effect indices of each input. + This approach, widely described in the paper by [saltelli1999]_, relies upon a + Fourier decomposition of the model response. Its key idea is to recast this + representation as a function of a *scalar* parameter :math:`s`, by defining + parametric curves :math:`s \mapsto x_i(s), i=1, \dots, n_X` exploring the + support of the input random vector :math:`\vect{X}`. + + Then the Fourier expansion of the model response is: + + .. math:: + + f(s) = \sum_{k \in \Zset^N} A_k cos(ks) + B_k sin(ks) + + where :math:`A_k` and :math:`B_k` are Fourier coefficients whose estimates are: + + .. math:: + + \hat{A}_k &= \frac{1}{N} \sum_{j=1}^N f(x_j^1,\dots,x_j^{N_X}) cos\left(\frac{2k\pi (j-1)}{N} \right) \quad , \quad -\frac{N}{2} \leq k \leq \frac{N}{2} \\ + \hat{B}_k &= \frac{1}{N} \sum_{j=1}^N f(x_j^1,\dots,x_j^{N_X}) sin\left(\frac{2k\pi (j-1)}{N} \right) \quad , \quad -\frac{N}{2} \leq k \leq \frac{N}{2} + + + The first order indices are estimated by: + + .. math:: + + \hat{S}_i = \frac{\hat{D}_i}{\hat{D}} + = \frac{\sum_{p=1}^M(\hat{A}_{p\omega_i}^2 + \hat{B}_{p\omega_i}^2)^2} + {\sum_{n=1}^{(N-1)/2}(\hat{A}_n^2 + \hat{B}_n^2)^2} + + and the total order indices by: + + .. math:: + + \hat{T}_i = 1 - \frac{\hat{D}_{-i}}{\hat{D}} + = 1 - \frac{\sum_{k=1}^{\omega_i/2}(\hat{A}_k^2 + \hat{B}_k^2)^2} + {\sum_{n=1}^{(N-1)/2}(\hat{A}_n^2 + \hat{B}_n^2)^2} + + where :math:`\hat{D}` is the total variance, :math:`\hat{D}_i` the portion + of :math:`D` arising from the uncertainty of the :math:`i^{th}` input and + :math:`\hat{D}_{-i}` is the part of the variance due to all the inputs + except the :math:`i^{th}` input. + + :math:`N` is the size of the sample using to compute the Fourier series and + :math:`M` is the interference factor. *Saltelli et al.* (1999) recommended to + set :math:`M` to a value in the range :math:`[4, 6]`. + :math:`\{\omega_i\}, \forall i=1, \dots, n_X` is a set of integer frequencies + assigned to each input :math:`X^i`. The frequency associated with the input + for which the sensitivity indices are computed, is set to the maximum admissible + frequency satisfying the Nyquist criterion (which ensures to avoid aliasing effects): + + .. math:: + + \omega_i = \frac{N - 1}{2M} + + In the paper by Saltelli et al. (1999), for high sample size, it is suggested + that :math:`16 \leq \omega_i/N_r \leq 64`. + +The computation relies on +`OpenTURNS capabilities `_. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass +from dataclasses import field +from typing import TYPE_CHECKING +from typing import Any +from typing import Final + +import matplotlib.pyplot as plt +import openturns as ot +from matplotlib.transforms import Affine2D +from numpy import array +from openturns import FFT +from strenum import PascalCaseStrEnum +from strenum import StrEnum + +from gemseo.core.mdo_functions.discipline_adapter import DisciplineAdapter +from gemseo.scenarios.doe_scenario import DOEScenario +from gemseo.uncertainty.sensitivity.base_sensitivity_analysis import ( + BaseSensitivityAnalysis, +) +from gemseo.uncertainty.sensitivity.base_sensitivity_analysis import ( + FirstOrderIndicesType, +) +from gemseo.utils.matplotlib_figure import save_show_figure_from_file_path_manager +from gemseo.utils.string_tools import filter_names +from gemseo.utils.string_tools import repr_variable + +if TYPE_CHECKING: + from collections.abc import Collection + from collections.abc import Iterable + from collections.abc import Sequence + from pathlib import Path + + from matplotlib.figure import Figure + + from gemseo.algos.parameter_space import ParameterSpace + from gemseo.core.discipline import Discipline + from gemseo.utils.string_tools import VariableType + +LOGGER = logging.getLogger(__name__) + + +class FAST(BaseSensitivityAnalysis): + """Sensitivity analysis based on the FAST (Fourier Amplitude Sensitivity Test) indices. + + Examples: + 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.sobol_analysis import SobolAnalysis + + configure_logger() + + discipline = IshigamiDiscipline() + uncertain_space = IshigamiSpace() + fast_analysis = FAST(disciplines=[discipline], parameter_space=uncertain_space, n_samples=400, output_names="y") + print(fast_analysis.indices) + """ + + @dataclass(frozen=True) + class SensitivityIndices: # noqa: D106 + first: FirstOrderIndicesType = field(default_factory=dict) + """The first-order Sobol' indices.""" + + total: FirstOrderIndicesType = field(default_factory=dict) + """The total order Sobol' indices.""" + + _indices: SensitivityIndices + + class Algorithm(PascalCaseStrEnum): + """The algorithms to estimate the Sobol' indices.""" + + FFT = "FFT" + + __ALGO_NAME_TO_CLASS: Final[dict[Algorithm, type]] = { + Algorithm.FFT: FFT, + } + """The mapping from the sensitivity algorithms to the OT classes.""" + + class Method(StrEnum): + """The names of the sensitivity methods.""" + + FIRST = "first" + """The first-order FAST index.""" + + TOTAL = "total" + """The total-order FAST index.""" + + def __init__( + self, + disciplines: Collection[Discipline], + parameter_space: ParameterSpace, + n_samples: int, + output_names: Iterable[str] = (), + formulation_name: str = "MDF", + **formulation_settings: Any, + ) -> None: + """Args: + disciplines: The discipline or disciplines to use for the analysis. + + parameter_space: The uncertain space containing the probability distribution of input parameters. + + n_samples: Number of samples to compute the FAST indices. + + output_names: The disciplines' outputs to be considered for the analysis. If empty, use all the outputs. + + formulation_name: The name of the :class:`.BaseMDOFormulation` to sample the disciplines. + + **formulation_settings: The settings of the :class:`.BaseMDOFormulation`. + + """ + self._indices = self.compute_indices( + disciplines, + parameter_space, + n_samples, + output_names, + formulation_name, + **formulation_settings, + ) + + def compute_indices( + self, + disciplines: Collection[Discipline], + parameter_space: ParameterSpace, + n_samples: int, + output_names: Iterable[str] = (), + formulation_name: str = "MDF", + **formulation_settings: Any, + ) -> SensitivityIndices: + scenario = self._create_scenario( + disciplines, + output_names, + formulation_name, + parameter_space, + **formulation_settings, + ) + disc = scenario.formulation.get_top_level_disciplines()[0] + disc.set_cache(disc.CacheType.MEMORY_FULL) + adapter = DisciplineAdapter( + parameter_space.variable_names, output_names, disc.default_input_data, disc + ) + + def evaluate(x_vect): + return [array(adapter.evaluate(array(x_vect))).tolist()] + + function = ot.PythonFunction( + parameter_space.dimension, len(output_names), evaluate + ) + self._first_order_indices = ot.FAST( + function, parameter_space.distribution.distribution, n_samples + ).getFirstOrderIndices() + self._total_order_indices = ot.FAST( + function, parameter_space.distribution.distribution, n_samples + ).getTotalOrderIndices() + self._all_indices = self.SensitivityIndices( + first=self._first_order_indices, + total=self._total_order_indices, + ) + return self._all_indices + + def _create_scenario( + self, + disciplines: Iterable[Discipline], + observable_names: Sequence[str], + formulation: str, + parameter_space: ParameterSpace, + ) -> DOEScenario: + """Create a DOE scenario to sample the disciplines. + + Args: + disciplines: The disciplines to sample. + observable_names: The names of the observables. + formulation: The name of the :class:`.MDOFormulation` + to sample the disciplines. + formulation_options: The options of the :class:`.MDOFormulation`. + parameter_space: A parameter space. + + Returns: + The DOE scenario to be used to sample the disciplines. + """ + scenario = DOEScenario( + disciplines, + observable_names, + parameter_space, + name=f"{self.__class__.__name__}SamplingPhase", + formulation_name=formulation, + ) + for _discipline in disciplines: + for output_name in observable_names: + if output_name in observable_names[1:]: + scenario.add_observable(output_name) + return scenario + + def plot( + self, + output: VariableType, + input_names: Iterable[str] = (), + title: str = "", + save: bool = True, + show: bool = False, + file_path: str | Path = "", + directory_path: str | Path = "", + file_name: str = "", + file_format: str = "", + sort: bool = True, + sort_by_total: bool = True, + ) -> Figure: + r"""Plot the first- and total-order Sobol' indices. + + For the :math:`i`-th input variable, + plot its first-order Sobol' index :math:`S_i^{1}` + and its total-order Sobol' index :math:`S_i^{T}` with dots + and their confidence intervals with vertical lines. + + The subtitle displays the standard deviation (StD) and the variance (Var) + of the output of interest. + + Args: + directory_path: The path to the directory where to save the plots. + file_name: The name of the file. + title: The title of the plot. + If empty, use a default one. + sort: Whether to sort the input variables by decreasing order. + sort_by_total: Whether to sort according to the total-order Sobol' indices + when ``sort`` is ``True``. + Otherwise, use the first-order Sobol' indices. + """ # noqa: D415 D417 + if not isinstance(output, tuple): + output = (output, 0) + + fig, ax = plt.subplots() + + indices = self.indices.total if sort_by_total else self.indices.first + output_name, output_component = output + indices = indices[output_name][output_component] + if sort: + names = [ + name + for name, _ in sorted( + indices.items(), key=lambda item: item[1].sum(), reverse=True + ) + ] + else: + names = indices.keys() + + names = filter_names(names, input_names) + + first_order_indices = self.indices.first[output_name][output_component] + total_order_indices = self.indices.total[output_name][output_component] + names_to_sizes = { + name: value.size for name, value in first_order_indices.items() + } + values_first_order = [ + first_order_indices[name][index] + for name in names + for index in range(names_to_sizes[name]) + ] + values_total_order = [ + total_order_indices[name][index] + for name in names + for index in range(names_to_sizes[name]) + ] + x_labels = [] + for name in names: + if names_to_sizes[name] == 1: + x_labels.append(name) + else: + size = names_to_sizes[name] + x_labels.extend([ + repr_variable(name, index, size) for index in range(size) + ]) + pretty_output_name = repr_variable( + output_name, + output_component, + len(self.indices.total[output_name]), + ) + if not title: + title = f"Sobol' indices for the output {pretty_output_name}" + variance = self.output_variances[output_name][output_component] + ax.set_title(f"{title}\nVar={variance:.1e} StD={variance**0.5:.1e}") + ax.set_axisbelow(True) + ax.grid() + + intervals = self.get_intervals() + intervals = intervals[output_name][output_component] + errorbar_options = {"marker": "o", "linestyle": "", "markersize": 7} + trans1 = Affine2D().translate(-0.01, 0.0) + ax.transData + trans2 = Affine2D().translate(+0.01, 0.0) + ax.transData + yerr = array([ + [ + first_order_indices[name][index] - intervals[name][0, index], + intervals[name][1, index] - first_order_indices[name][index], + ] + for name in names + for index in range(names_to_sizes[name]) + ]).T + + ax.errorbar( + x_labels, + values_first_order, + yerr=yerr, + label="First order", + transform=trans2, + **errorbar_options, + ) + intervals = self.get_intervals(False) + intervals = intervals[output_name][output_component] + yerr = array([ + [ + total_order_indices[name][index] - intervals[name][0, index], + intervals[name][1, index] - total_order_indices[name][index], + ] + for name in names + for index in range(names_to_sizes[name]) + ]).T + ax.errorbar( + x_labels, + values_total_order, + yerr, + label="Total order", + transform=trans1, + **errorbar_options, + ) + ax.legend(loc="lower left") + save_show_figure_from_file_path_manager( + fig, + self._file_path_manager if save else None, + show=show, + file_path=file_path, + file_name=file_name, + file_format=file_format, + directory_path=directory_path, + ) + return fig diff --git a/tests/uncertainty/sensitivity/test_fast.py b/tests/uncertainty/sensitivity/test_fast.py new file mode 100644 index 0000000000000000000000000000000000000000..622f1b3650682dd92d7f90e13a49f2b254b6d7c7 --- /dev/null +++ b/tests/uncertainty/sensitivity/test_fast.py @@ -0,0 +1,56 @@ +# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License version 3 as published by the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com +# +# This work is licensed under a BSD 0-Clause License. +# +# Permission to use, copy, modify, and/or distribute this software +# for any purpose with or without fee is hereby granted. +# +# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL +# WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL +# THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, +# OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING +# FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, +# NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION +# WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. +# Contributors: +# INITIAL AUTHORS - initial API and implementation and/or initial +# documentation +# :author: Gokulnath Ramachandran +# OTHER AUTHORS - MACROSCOPIC CHANGES +""" +FAST analysis +=============== +""" + +from __future__ import annotations + +from gemseo.problems.uncertainty.ishigami.ishigami_discipline import IshigamiDiscipline +from gemseo.problems.uncertainty.ishigami.ishigami_space import IshigamiSpace +from gemseo.uncertainty.sensitivity.fast_analysis import FAST + + +def test_fast_method(): + discipline = IshigamiDiscipline() + uncertain_space = IshigamiSpace() + FAST( + disciplines=[discipline], + parameter_space=uncertain_space, + n_samples=400, + output_names="y", + )