From 80f360958368b44cbefe9b2e57f312b15d3a0aee Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 26 Jul 2023 16:02:26 -0400 Subject: [PATCH 1/9] added test for seed_u is None in FORM.run() --- tests/unit_tests/reliability/test_form.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/unit_tests/reliability/test_form.py b/tests/unit_tests/reliability/test_form.py index dadbd457..6d376cb7 100644 --- a/tests/unit_tests/reliability/test_form.py +++ b/tests/unit_tests/reliability/test_form.py @@ -43,6 +43,15 @@ def test_seeds_x_is_none(setup): np.testing.assert_allclose(form_obj.failure_probability, 0.0126, rtol=1e-02) +def test_seed_u_is_none(setup): + """ToDo: Fix FORM.run(seed_x=numpy_array) to pass this test""" + distributions = [Normal(loc=200, scale=20), Normal(loc=150, scale=10)] + form = FORM(distributions=distributions, runmodel_object=setup) + seed_x = np.array([225, 140]) + form.run(seed_x=seed_x) + np.testing.assert_allclose(form.failure_probability, 0.0126, rtol=1e-02) + + def test_tol1_is_not_none(setup): path = os.path.abspath(os.path.dirname(__file__)) os.chdir(path) From ec3f6e9949a9f3c3bf2b50ca6957afb62c348b14 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Thu, 27 Jul 2023 09:21:14 -0400 Subject: [PATCH 2/9] passes test for seed_u_is_none in FORM.run() --- src/UQpy/reliability/taylor_series/FORM.py | 2 +- tests/unit_tests/reliability/test_form.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index c10c895c..acd94e10 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -156,7 +156,7 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda elif seed_u is None and seed_x is not None: self.nataf_object.run(samples_x=seed_x.reshape(1, -1), jacobian=False) seed_z = self.nataf_object.samples_z - seed = Decorrelate(seed_z, self.nataf_object.corr_z) + seed = Decorrelate(seed_z, self.nataf_object.corr_z).samples_u elif seed_u is not None and seed_x is None: seed = np.squeeze(seed_u) else: diff --git a/tests/unit_tests/reliability/test_form.py b/tests/unit_tests/reliability/test_form.py index 6d376cb7..266acf87 100644 --- a/tests/unit_tests/reliability/test_form.py +++ b/tests/unit_tests/reliability/test_form.py @@ -44,7 +44,6 @@ def test_seeds_x_is_none(setup): def test_seed_u_is_none(setup): - """ToDo: Fix FORM.run(seed_x=numpy_array) to pass this test""" distributions = [Normal(loc=200, scale=20), Normal(loc=150, scale=10)] form = FORM(distributions=distributions, runmodel_object=setup) seed_x = np.array([225, 140]) From ee0fa3e1ab736f163e63c10c7724cf6971e3bd16 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Thu, 27 Jul 2023 09:30:02 -0400 Subject: [PATCH 3/9] added Connor to development team --- README.rst | 4 +--- docs/source/index.rst | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index 919d8759..6f9d6916 100644 --- a/README.rst +++ b/README.rst @@ -6,8 +6,6 @@ .. |AzureDevops| image:: https://img.shields.io/azure-devops/build/UQpy/5ce1851f-e51f-4e18-9eca-91c3ad9f9900/1?style=plastic :alt: Azure DevOps builds .. |PyPIdownloads| image:: https://img.shields.io/pypi/dm/UQpy?style=plastic :alt: PyPI - Downloads .. |PyPI| image:: https://img.shields.io/pypi/v/UQpy?style=plastic :alt: PyPI -.. |Binder| image:: https://mybinder.org/badge_logo.svg - :target: https://mybinder.org/v2/gh/SURGroup/UQpy/master .. |bear-ified| image:: https://raw.githubusercontent.com/beartype/beartype-assets/main/badge/bear-ified.svg :align: top @@ -32,7 +30,7 @@ Uncertainty Quantification with python (UQpy) + + + | | Ketson RM dos Santos, Katiana Kontolati, Dimitris Loukrezis, | + + + -| | Promit Chakroborty, Lukáš Novák, Andrew Solanto | +| | Promit Chakroborty, Lukáš Novák, Andrew Solanto, Connor Krill | +-----------------------+------------------------------------------------------------------+ | **Contributors:** | Michael Gardner, Prateek Bhustali, Julius Schultz, Ulrich Römer | +-----------------------+------------------------------------------------------------------+ diff --git a/docs/source/index.rst b/docs/source/index.rst index c3a61ae6..fc7af386 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -18,7 +18,7 @@ as a set of modules centered around core capabilities in Uncertainty Quantificat + + + | | Ketson RM dos Santos, Katiana Kontolati, Dimitris Loukrezis, | + + + -| | Promit Chakroborty, Lukáš Novák, Andrew Solanto | +| | Promit Chakroborty, Lukáš Novák, Andrew Solanto, Connor Krill | +-----------------------+------------------------------------------------------------------+ | **Contributors:** | Michael Gardner, Prateek Bhustali, Julius Schultz, Ulrich Römer | +-----------------------+------------------------------------------------------------------+ From 174d578c17c4811415e47b4d3594bc6addf574ab Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Thu, 27 Jul 2023 17:13:43 +0300 Subject: [PATCH 4/9] Minor changes to fix reliability examples --- ...nction_3d.py => FORM_linear_function_3d.py} | 14 +++++++------- docs/code/reliability/sorm/local_model4.py | 8 ++++++++ docs/code/reliability/sorm/local_pfn.py | 18 ++++++++---------- .../sorm/plot_SORM_nonlinear_function.py | 6 +++--- 4 files changed, 26 insertions(+), 20 deletions(-) rename docs/code/reliability/form/{plot_FORM_linear_function_3d.py => FORM_linear_function_3d.py} (71%) create mode 100644 docs/code/reliability/sorm/local_model4.py diff --git a/docs/code/reliability/form/plot_FORM_linear_function_3d.py b/docs/code/reliability/form/FORM_linear_function_3d.py similarity index 71% rename from docs/code/reliability/form/plot_FORM_linear_function_3d.py rename to docs/code/reliability/form/FORM_linear_function_3d.py index 055d7272..2e492e0f 100644 --- a/docs/code/reliability/form/plot_FORM_linear_function_3d.py +++ b/docs/code/reliability/form/FORM_linear_function_3d.py @@ -35,13 +35,13 @@ dist3 = Normal(loc=4., scale=0.4) model = PythonModel(model_script='local_pfn.py', model_object_name="example3",) -RunModelObject3 = RunModel(model=model) +run_model = RunModel(model=model) -Z0 = FORM(distributions=[dist1, dist2, dist3], runmodel_object=RunModelObject3) -Z0.run() +form = FORM(distributions=[dist1, dist2, dist3], runmodel_object=run_model) +form.run() -print('Design point in standard normal space: %s' % Z0.design_point_u) -print('Design point in original space: %s' % Z0.design_point_x) -print('Hasofer-Lind reliability index: %s' % Z0.beta) -print('FORM probability of failure: %s' % Z0.failure_probability) +print('Design point in standard normal space: %s' % form.design_point_u) +print('Design point in original space: %s' % form.design_point_x) +print('Hasofer-Lind reliability index: %s' % form.beta) +print('FORM probability of failure: %s' % form.failure_probability) diff --git a/docs/code/reliability/sorm/local_model4.py b/docs/code/reliability/sorm/local_model4.py new file mode 100644 index 00000000..6c0808cc --- /dev/null +++ b/docs/code/reliability/sorm/local_model4.py @@ -0,0 +1,8 @@ +import numpy as np + + +def example4(samples=None): + g = np.zeros(samples.shape[0]) + for i in range(samples.shape[0]): + g[i] = samples[i, 0] * samples[i, 1] - 80 + return g \ No newline at end of file diff --git a/docs/code/reliability/sorm/local_pfn.py b/docs/code/reliability/sorm/local_pfn.py index db3584d0..2903bf25 100644 --- a/docs/code/reliability/sorm/local_pfn.py +++ b/docs/code/reliability/sorm/local_pfn.py @@ -6,14 +6,15 @@ """ import numpy as np + def example1(samples=None): g = np.zeros(samples.shape[0]) - for i in range(samples.shape[0]): + for i in range(samples.shape[0]): R = samples[i, 0] S = samples[i, 1] g[i] = R - S return g - + def example2(samples=None): import numpy as np @@ -21,18 +22,15 @@ def example2(samples=None): beta = 3.0902 g = np.zeros(samples.shape[0]) for i in range(samples.shape[0]): - g[i] = -1/np.sqrt(d) * (samples[i, 0] + samples[i, 1]) + beta + g[i] = -1 / np.sqrt(d) * (samples[i, 0] + samples[i, 1]) + beta return g def example3(samples=None): g = np.zeros(samples.shape[0]) for i in range(samples.shape[0]): - g[i] = 6.2*samples[i, 0] - samples[i, 1]*samples[i, 2]**2 + g[i] = 6.2 * samples[i, 0] - samples[i, 1] * samples[i, 2] ** 2 return g - -def example4(samples=None): - g = np.zeros(samples.shape[0]) - for i in range(samples.shape[0]): - g[i] = samples[i, 0]*samples[i, 1] - 80 - return g \ No newline at end of file + + + diff --git a/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py b/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py index 25538613..05b2f32d 100644 --- a/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py +++ b/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py @@ -34,13 +34,13 @@ dist1 = Normal(loc=20., scale=2) dist2 = Lognormal(s=s, loc=0.0, scale=scale) -model = PythonModel(model_script='local_pfn.py', model_object_name="example4",) +model = PythonModel(model_script='local_model4.py', model_object_name="example4") RunModelObject4 = RunModel(model=model) form = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject4) form.run() -Q0 = SORM(form_object=form) +sorm = SORM(form_object=form) # print results -print('SORM probability of failure: %s' % Q0.failure_probability) +print('SORM probability of failure: %s' % sorm.failure_probability) From 7a5bd248576dcced9b9be77f0c003d7564ba7d6b Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Thu, 27 Jul 2023 17:47:44 +0300 Subject: [PATCH 5/9] Corrects equation display in plot_pce_sparsity_lars.py --- .../surrogates/pce/plot_pce_sparsity_lars.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/code/surrogates/pce/plot_pce_sparsity_lars.py b/docs/code/surrogates/pce/plot_pce_sparsity_lars.py index 864e9110..d3bf7363 100644 --- a/docs/code/surrogates/pce/plot_pce_sparsity_lars.py +++ b/docs/code/surrogates/pce/plot_pce_sparsity_lars.py @@ -15,8 +15,7 @@ # %% import numpy as np -import math -import numpy as np + from UQpy.distributions import Uniform, JointIndependent from UQpy.surrogates import * @@ -24,7 +23,8 @@ # %% md # # We then define the Ishigami function, which reads: -# :math:` f(x_1, x_2, x_3) = \sin(x_1) + a \sin^2(x_2) + b x_3^4 \sin(x_1)` +# +# .. math:: f(x_1, x_2, x_3) = \sin(x_1) + a \sin^2(x_2) + b x_3^4 \sin(x_1) # %% @@ -41,7 +41,7 @@ def ishigami(xx): # %% md # -# The Ishigami function has three indepdent random inputs, which are uniformly distributed in +# The Ishigami function has three independent random inputs, which are uniformly distributed in # interval :math:`[-\pi, \pi]`. # %% @@ -70,7 +70,7 @@ def ishigami(xx): # # where :math:`N` is the number of random inputs (here, :math:`N=3`). # -# Note that the size of the basis is highly dependent both on :math:`N` and :math:`P:math:`. It is generally advisable +# Note that the size of the basis is highly dependent both on :math:`N` and :math:`P`. It is generally advisable # that the experimental design has :math:`2-10` times more data points than the number of PCE polynomials. This might # lead to curse of dimensionality and thus we will utilize the best model selection algorithm based on # Least Angle Regression. @@ -102,8 +102,8 @@ def ishigami(xx): # %% md # -# We now fit the PCE coefficients by solving a regression problem. Here we opt for the _np.linalg.lstsq_ method, -# which is based on the _dgelsd_ solver of LAPACK. This original PCE class will be used for further selection of +# We now fit the PCE coefficients by solving a regression problem. Here we opt for the :code:`_np.linalg.lstsq_` method, +# which is based on the :code:`_dgelsd_` solver of LAPACK. This original PCE class will be used for further selection of # the best basis functions. # %% @@ -238,7 +238,7 @@ def ishigami(xx): # %% md # # In case of high-dimensional input and/or high :math:P` it is also beneficial to reduce the TD basis set by hyperbolic -# trunction. The hyperbolic truncation reduces higher-order interaction terms in dependence to parameter :math:`q` in +# truncation. The hyperbolic truncation reduces higher-order interaction terms in dependence to parameter :math:`q` in # interval :math:`(0,1)`. The set of multi indices :math:`\alpha` is reduced as follows: # # :math:`\alpha\in \mathbb{N}^{N}: || \boldsymbol{\alpha}||_q \equiv \Big( \sum_{i=1}^{N} \alpha_i^q \Big)^{1/q} \leq P` From 81335384ef26af127077df38e9aeb1fb407ce9dc Mon Sep 17 00:00:00 2001 From: connor-krill Date: Thu, 17 Aug 2023 13:22:30 -0400 Subject: [PATCH 6/9] grammatical and notation edits to make documentation more consistent with literature references --- docs/source/reliability/index.rst | 15 ++++++++------- docs/source/reliability/subset.rst | 24 ++++++++++++------------ 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/source/reliability/index.rst b/docs/source/reliability/index.rst index d3010039..51362cfb 100644 --- a/docs/source/reliability/index.rst +++ b/docs/source/reliability/index.rst @@ -2,19 +2,20 @@ Reliability =========== Reliability of a system refers to the assessment of its probability of failure (i.e the system no longer satisfies some -performance measures), given the model uncertainty in the structural, environmental and load parameters. Given a vector -of random variables :math:`\textbf{X}=\{X_1, X_2, \ldots, X_n\} \in \mathcal{D}_\textbf{X}\subset \mathbb{R}^n`, where -:math:`\mathcal{D}` is the domain of interest and :math:`f_{\textbf{X}}(\textbf{x})` is its joint probability density -function then, the probability that the system will fail is defined as +performance measure), given the model uncertainty in the structural, environmental and load parameters. Given a vector +of random variables :math:`\textbf{X}=[X_1, X_2, \ldots, X_n]^T \in \mathcal{D}_\textbf{X}\subset \mathbb{R}^n`, where +:math:`\mathcal{D}_\textbf{X}` is the domain of interest and :math:`f_{\textbf{X}}(\textbf{x})` is its joint probability density +function, then the probability that the system will fail is defined as -.. math:: P_f =\mathbb{P}(g(\textbf{X}) \leq 0) = \int_{D_f} f_{\textbf{X}}(\textbf{x})d\textbf{x} = \int_{\{\textbf{X}:g(\textbf{X})\leq 0 \}} f_{\textbf{X}}(\textbf{x})d\textbf{x} +.. math:: P_f =\mathbb{P}(g(\textbf{X}) \leq 0) = \int_{\mathcal{D}_f} f_{\textbf{X}}(\textbf{x})d\textbf{x} = \int_{\{\textbf{X}:g(\textbf{X})\leq 0 \}} f_{\textbf{X}}(\textbf{x})d\textbf{x} -where :math:`g(\textbf{X})` is the so-called performance function. The reliability problem is often formulated in the +where :math:`g(\textbf{X})` is the so-called performance function and :math:`\mathcal{D}_f=\{\textbf{X}:g(\textbf{X})\leq 0 \}` is the failure domain. +The reliability problem is often formulated in the standard normal space :math:`\textbf{U}\sim \mathcal{N}(\textbf{0}, \textbf{I}_n)`, which means that a nonlinear isoprobabilistic transformation from the generally non-normal parameter space -:math:`\textbf{X}\sim f_{\textbf{X}}(\cdot)` to the standard normal is required (see the :py:mod:`.transformations` module). +:math:`\textbf{X}\sim f_{\textbf{X}}(\cdot)` to the standard normal space is required (see the :py:mod:`.transformations` module). The performance function in the standard normal space is denoted :math:`G(\textbf{U})`. :py:mod:`.UQpy` does not require this transformation and can support reliability analysis for problems with arbitrarily distributed parameters. diff --git a/docs/source/reliability/subset.rst b/docs/source/reliability/subset.rst index 5da7def8..555b76af 100644 --- a/docs/source/reliability/subset.rst +++ b/docs/source/reliability/subset.rst @@ -2,27 +2,27 @@ Subset Simulation ------------------- In the subset simulation method :cite:`SubsetSimulation` the probability of failure :math:`P_f` is approximated by a product of probabilities -of more frequent events. That is, the failure event :math:`G = \{\textbf{x} \in \mathbb{R}^n:G(\textbf{x}) \leq 0\}`, +of more frequent events. That is, the failure event :math:`G = \{\textbf{X} \in \mathbb{R}^n:g(\textbf{X}) \leq 0\}`, is expressed as the of union of `M` nested intermediate events :math:`G_1,G_2,\cdots,G_M` such that :math:`G_1 \supset G_2 \supset \cdots \supset G_M`, and :math:`G = \cap_{i=1}^{M} G_i`. The intermediate failure events -are defined as :math:`G_i=\{G(\textbf{x})\le b_i\}`, where :math:`b_1>b_2>\cdots>b_i=0` are positive thresholds selected -such that each conditional probability :math:`P(G_i | G_{i-1}), ~i=2,3,\cdots,M-1` equals a target probability value +are defined as :math:`G_i=\{g(\textbf{X})\le b_i\}`, where :math:`b_1>b_2>\cdots>b_M=0` are non-negative thresholds selected +such that each conditional probability :math:`P(G_{i+1} | G_{i}),\ i=1,2,\cdots,M-1` equals a target probability value :math:`p_0`. The probability of failure :math:`P_f` is estimated as: -.. math:: P_f = P\left(\cap_{i=1}^M G_i\right) = P(G_1)\prod_{i=2}^M P(G_i | G_{i-1}) +.. math:: P_f = P\left(\bigcap_{i=1}^M G_i\right) = P(G_1)\prod_{i=1}^{M-1} P(G_{i+1} | G_{i}) where the probability :math:`P(G_1)` is computed through Monte Carlo simulations. In order to estimate the conditional -probabilities :math:`P(G_i|G_{i-1}),~j=2,3,\cdots,M` generation of Markov Chain Monte Carlo (MCMC) samples from the -conditional pdf :math:`p_{\textbf{U}}(\textbf{u}|G_{i-1})` is required. In the context of subset simulation, the Markov +probabilities :math:`P(G_{i+1}|G_i),~i=1,2,\cdots,M-1` generation of Markov Chain Monte Carlo (MCMC) samples from the +conditional pdf :math:`p_{\textbf{X}}(\textbf{x}|G_i)` is required. In the context of subset simulation, the Markov chains are constructed through a two-step acceptance/rejection criterion. Starting from a Markov chain state -:math:`\textbf{x}` and a proposal distribution :math:`q(\cdot|\textbf{x})`, a candidate sample :math:`\textbf{w}` is -generated. In the first stage, the sample :math:`\textbf{w}` is accepted/rejected with probability +:math:`\textbf{X}` and a proposal distribution :math:`q(\cdot|\textbf{X})`, a candidate sample :math:`\textbf{W}` is +generated. In the first stage, the sample :math:`\textbf{W}` is accepted/rejected with probability -.. math:: \alpha=\min\bigg\{1, \frac{p(\textbf{w})q(\textbf{x}|\textbf{w})}{p(\textbf{x})q(\textbf{w}|\textbf{x})}\bigg\} +.. math:: \alpha=\min\bigg\{1, \frac{p_\textbf{X}(\textbf{w})q(\textbf{x}|\textbf{W})}{p_\textbf{X}(\textbf{x})q(\textbf{w}|\textbf{X})}\bigg\} -and in the second stage is accepted/rejected based on whether the sample belongs to the failure region :math:`G_j`. -:class:`.SubSetSimulation` can be used with any of the available (or custom) :class:`.MCMC` classes in the -:py:mod:`sampling` module. +and in the second stage is accepted/rejected based on whether the sample belongs to the failure region :math:`G_i`. +:class:`.SubsetSimulation` can be used with any of the available (or custom) :class:`.MCMC` classes in the +:py:mod:`Sampling` module. SubsetSimulation Class ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From a75aac526ec918fe4a95b7df38e059cc99acb63b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 18 Aug 2023 13:07:52 -0400 Subject: [PATCH 7/9] cleaned up notation and attribute documentation --- docs/source/reliability/subset.rst | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/docs/source/reliability/subset.rst b/docs/source/reliability/subset.rst index 555b76af..be50bfc9 100644 --- a/docs/source/reliability/subset.rst +++ b/docs/source/reliability/subset.rst @@ -24,26 +24,24 @@ and in the second stage is accepted/rejected based on whether the sample belongs :class:`.SubsetSimulation` can be used with any of the available (or custom) :class:`.MCMC` classes in the :py:mod:`Sampling` module. -SubsetSimulation Class -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - The :class:`.SubsetSimulation` class is imported using the following command: >>> from UQpy.reliability.SubsetSimulation import SubsetSimulation Methods """"""" + .. autoclass:: UQpy.reliability.SubsetSimulation Attributes """""""""" -.. autoattribute:: UQpy.reliability.SubsetSimulation.samples -.. autoattribute:: UQpy.reliability.SubsetSimulation.performance_function_per_level -.. autoattribute:: UQpy.reliability.SubsetSimulation.performance_threshold_per_level + +.. autoattribute:: UQpy.reliability.SubsetSimulation.dependent_chains_CoV .. autoattribute:: UQpy.reliability.SubsetSimulation.failure_probability .. autoattribute:: UQpy.reliability.SubsetSimulation.independent_chains_CoV -.. autoattribute:: UQpy.reliability.SubsetSimulation.dependent_chains_CoV - +.. autoattribute:: UQpy.reliability.SubsetSimulation.performance_function_per_level +.. autoattribute:: UQpy.reliability.SubsetSimulation.performance_threshold_per_level +.. autoattribute:: UQpy.reliability.SubsetSimulation.samples Examples """""""""" From adad1d0cb1ee73b51bc70b8e8bf3c8811ff2f3cf Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 18 Aug 2023 13:08:06 -0400 Subject: [PATCH 8/9] deleted unused file --- docs/code/reliability/subset_simulation/pfn.py | 11 ----------- 1 file changed, 11 deletions(-) delete mode 100644 docs/code/reliability/subset_simulation/pfn.py diff --git a/docs/code/reliability/subset_simulation/pfn.py b/docs/code/reliability/subset_simulation/pfn.py deleted file mode 100644 index 9dbd2179..00000000 --- a/docs/code/reliability/subset_simulation/pfn.py +++ /dev/null @@ -1,11 +0,0 @@ -""" - -Auxiliary File -====================================================================== - -""" -import numpy as np - - -def run_python_model(samples, b_eff, d): - return [b_eff * np.sqrt(d) - np.sum(samples[i, :]) for i in range(samples.shape[0])] From 6359c16806620eabbb1664ee193d5eec274450bf Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 18 Aug 2023 13:09:53 -0400 Subject: [PATCH 9/9] cleaned up and added docstrings to public and hidden methods, renamed internal variables for clarity and consistency with literature --- src/UQpy/reliability/SubsetSimulation.py | 240 ++++++++++++----------- 1 file changed, 124 insertions(+), 116 deletions(-) diff --git a/src/UQpy/reliability/SubsetSimulation.py b/src/UQpy/reliability/SubsetSimulation.py index c0d5f663..99e8989d 100644 --- a/src/UQpy/reliability/SubsetSimulation.py +++ b/src/UQpy/reliability/SubsetSimulation.py @@ -1,9 +1,6 @@ -import copy -import logging - -from UQpy.run_model import RunModel -from UQpy.utilities.Utilities import process_random_state from UQpy.sampling import * +from UQpy.utilities.Utilities import process_random_state +from UQpy.utilities.ValidationTypes import PositiveInteger class SubsetSimulation: @@ -15,25 +12,25 @@ def __init__( sampling: MCMC, samples_init: np.ndarray = None, conditional_probability: Annotated[Union[float, int], Is[lambda x: 0 <= x <= 1]] = 0.1, - nsamples_per_subset: int = 1000, - max_level: int = 10, + nsamples_per_subset: PositiveInteger = 1000, + max_level: PositiveInteger = 10, ): """ Perform Subset Simulation to estimate probability of failure. This class estimates probability of failure for a user-defined model using Subset Simulation. The class can - use one of several mcmc algorithms to draw conditional samples. + use one of several MCMC algorithms to draw conditional samples. :param runmodel_object: The computational model. It should be of type :class:`.RunModel`. :param sampling: Specifies the :class:`MCMC` algorithm. :param samples_init: A set of samples from the specified probability distribution. These are the samples from the original distribution. They are not conditional samples. The samples must be an array of size :code:`samples_number_per_subset x dimension`. - If `samples_init` is not specified, the :class:`.SubsetSimulation` class will use the :class:`.MCMC` class - created with the aid of `sampling` to draw the initial samples. - :param conditional_probability: Conditional probability for each conditional level. - :param nsamples_per_subset: Number of samples to draw in each conditional level. - :param max_level: Maximum number of allowable conditional levels. + If :code:`samples_init` is not specified, the :class:`.SubsetSimulation` class will use the :class:`.MCMC` class + created with the aid of :code:`sampling` to draw the initial samples. + :param conditional_probability: Conditional probability for each conditional level. Default: :math:`0.1` + :param nsamples_per_subset: Number of samples to draw in each conditional level. Default: :math:`1000` + :param max_level: Maximum number of allowable conditional levels. Default: :math:`10` """ # Initialize other attributes self._sampling_class = sampling @@ -44,112 +41,105 @@ def __init__( self.nsamples_per_subset = nsamples_per_subset self.max_level = max_level self.logger = logging.getLogger(__name__) - self.mcmc_objects = [sampling] - self.samples: list = list() - """A list of arrays containing the samples in each conditional level. The size of the list is equal to the - number of levels.""" + self.dependent_chains_CoV: float = None + """Coefficient of variation of the probability of failure estimate with dependent chains.""" + self.failure_probability: float = None + """Probability of failure estimate.""" + self.independent_chains_CoV: float = None + """Coefficient of variation of the probability of failure estimate assuming independent chains.""" self.performance_function_per_level: list = [] """A list of arrays containing the evaluation of the performance function at each sample in each conditional level. The size of the list is equal to the number of levels.""" self.performance_threshold_per_level: list = [] """Threshold value of the performance function for each conditional level. The size of the list is equal to the number of levels.""" + self.samples: list = [] + """A list of arrays containing the samples in each conditional level. The size of the list is equal to the + number of levels.""" self.logger.info("UQpy: Running Subset Simulation with mcmc of type: " + str(type(sampling))) - self.failure_probability: float = None - """Probability of failure estimate.""" - self.independent_chains_CoV: float = None - """Coefficient of variation of the probability of failure estimate assuming independent chains.""" - self.dependent_chains_CoV: float = None - """Coefficient of variation of the probability of failure estimate with dependent chains.""" - [self.failure_probability, self.independent_chains_CoV, self.dependent_chains_CoV] = self._run() - self.logger.info("UQpy: Subset Simulation Complete!") - # The run function executes the chosen subset simulation algorithm def _run(self): - """ - Execute subset simulation + """Execute subset simulation This is an instance method that runs subset simulation. It is automatically called when the :class:`.SubsetSimulation` class is instantiated. + + :return: failure_probability, probability_cov_independent, probability_cov_dependent """ - step = 0 + conditional_level = 0 n_keep = int(self.conditional_probability * self.nsamples_per_subset) - d12 = [] - d22 = [] + independent_chain_cov_squared = [] + dependent_chain_cov_squared = [] # Generate the initial samples - Level 0 # Here we need to make sure that we have good initial samples from the target joint density. if self.samples_init is None: self.logger.warning( - "UQpy: You have not provided initial samples.\n Subset simulation is highly sensitive " - "to the initial sample set. It is recommended that the user either:\n" - "- Provide an initial set of samples (samples_init) known to follow the distribution; " - "or\n - Provide a robust mcmc object that will draw independent initial samples from " - "the distribution." + "UQpy: You have not provided initial samples." + "\n Subset simulation is highly sensitive to the initial sample set. It is recommended that the user either:" + "\n- Provide an initial set of samples (samples_init) known to follow the joint distribution of the parameters; or" + "\n- Provide a robust mcmc object that will draw independent initial samples from the joint distribution of the parameters." ) self.mcmc_objects[0].run(nsamples=self.nsamples_per_subset) self.samples.append(self.mcmc_objects[0].samples) else: self.samples.append(self.samples_init) - # Run the model for the initial samples, sort them by their performance function, and identify the - # conditional level - self.runmodel_object.run(samples=np.atleast_2d(self.samples[step])) + # Run the model with initial samples, sort by their performance function, and identify the conditional level + self.runmodel_object.run(samples=np.atleast_2d(self.samples[conditional_level])) self.performance_function_per_level.append(np.squeeze(self.runmodel_object.qoi_list)) - g_ind = np.argsort(self.performance_function_per_level[step]) - self.performance_threshold_per_level.append(self.performance_function_per_level[step][g_ind[n_keep - 1]]) + g_ind = np.argsort(self.performance_function_per_level[conditional_level]) + self.performance_threshold_per_level.append(self.performance_function_per_level[conditional_level][g_ind[n_keep - 1]]) # Estimate coefficient of variation of conditional probability of first level - d1, d2 = self._compute_coefficient_of_variation(step) - d12.append(d1 ** 2) - d22.append(d2 ** 2) + independent_intermediate_cov, dependent_intermediate_cov = self._compute_intermediate_cov(conditional_level) + independent_chain_cov_squared.append(independent_intermediate_cov ** 2) + dependent_chain_cov_squared.append(dependent_intermediate_cov ** 2) self.logger.info("UQpy: Subset Simulation, conditional level 0 complete.") - while self.performance_threshold_per_level[step] > 0 and step < self.max_level: - - # Increment the conditional level - step += 1 + while self.performance_threshold_per_level[conditional_level] > 0 and conditional_level < self.max_level: + conditional_level += 1 # Increment the conditional level # Initialize the samples and the performance function at the next conditional level - self.samples.append(np.zeros_like(self.samples[step - 1])) - self.samples[step][:n_keep] = self.samples[step - 1][g_ind[0:n_keep], :] - self.performance_function_per_level.append(np.zeros_like(self.performance_function_per_level[step - 1])) - self.performance_function_per_level[step][:n_keep] = self.performance_function_per_level[step - 1][g_ind[:n_keep]] + self.samples.append(np.zeros_like(self.samples[conditional_level - 1])) + self.samples[conditional_level][:n_keep] = self.samples[conditional_level - 1][g_ind[0:n_keep], :] + self.performance_function_per_level.append(np.zeros_like(self.performance_function_per_level[conditional_level - 1])) + self.performance_function_per_level[conditional_level][:n_keep] = self.performance_function_per_level[conditional_level - 1][g_ind[:n_keep]] # Unpack the attributes new_sampler = copy.deepcopy(self._sampling_class) - new_sampler.seed = np.atleast_2d(self.samples[step][:n_keep, :]) - new_sampler.n_chains = len(np.atleast_2d(self.samples[step][:n_keep, :])) + new_sampler.seed = np.atleast_2d(self.samples[conditional_level][:n_keep, :]) + new_sampler.n_chains = len(np.atleast_2d(self.samples[conditional_level][:n_keep, :])) new_sampler.random_state = process_random_state(self._random_state) self.mcmc_objects.append(new_sampler) # Set the number of samples to propagate each chain (n_prop) in the conditional level - n_prop_test = (self.nsamples_per_subset / self.mcmc_objects[step].n_chains) + n_prop_test = self.nsamples_per_subset / self.mcmc_objects[conditional_level].n_chains if n_prop_test.is_integer(): - n_prop = (self.nsamples_per_subset // self.mcmc_objects[step].n_chains) + n_prop = self.nsamples_per_subset // self.mcmc_objects[conditional_level].n_chains else: raise AttributeError( - "UQpy: The number of samples per subset (nsamples_per_ss) must be an integer multiple of " - "the number of mcmc chains.") + "UQpy: The number of samples per subset (nsamples_per_subset) must be an integer multiple of " + "the number of MCMC chains.") # Propagate each chain n_prop times and evaluate the model to accept or reject. for i in range(n_prop - 1): # Propagate each chain if i == 0: - self.mcmc_objects[step].run(nsamples=2 * self.mcmc_objects[step].n_chains) + self.mcmc_objects[conditional_level].run(nsamples=2 * self.mcmc_objects[conditional_level].n_chains) else: - self.mcmc_objects[step].run(nsamples=self.mcmc_objects[step].n_chains) + self.mcmc_objects[conditional_level].run(nsamples=self.mcmc_objects[conditional_level].n_chains) # Decide whether a new simulation is needed for each proposed state - a = self.mcmc_objects[step].samples[i * n_keep : (i + 1) * n_keep, :] - b = self.mcmc_objects[step].samples[(i + 1) * n_keep : (i + 2) * n_keep, :] + a = self.mcmc_objects[conditional_level].samples[i * n_keep : (i + 1) * n_keep, :] + b = self.mcmc_objects[conditional_level].samples[(i + 1) * n_keep : (i + 2) * n_keep, :] test1 = np.equal(a, b) test = np.logical_and(test1[:, 0], test1[:, 1]) @@ -159,66 +149,72 @@ def _run(self): ind_true = [i for i, val in enumerate(test) if val] # Do not run the model for those samples where the mcmc state remains unchanged. - self.samples[step][[x + (i + 1) * n_keep for x in ind_true], :] = \ - self.mcmc_objects[step].samples[ind_true, :] - self.performance_function_per_level[step][[x + (i + 1) * n_keep for x in ind_true]] = \ - self.performance_function_per_level[step][ind_true] + self.samples[conditional_level][[x + (i + 1) * n_keep for x in ind_true], :] = \ + self.mcmc_objects[conditional_level].samples[ind_true, :] + self.performance_function_per_level[conditional_level][[x + (i + 1) * n_keep for x in ind_true]] = \ + self.performance_function_per_level[conditional_level][ind_true] # Run the model at each of the new sample points - x_run = self.mcmc_objects[step].samples[[x + (i + 1) * n_keep for x in ind_false], :] + x_run = self.mcmc_objects[conditional_level].samples[[x + (i + 1) * n_keep for x in ind_false], :] if x_run.size != 0: self.runmodel_object.run(samples=x_run) # Temporarily save the latest model runs - g_temp = np.asarray(self.runmodel_object.qoi_list[-len(x_run) :]) + response_function_values = np.asarray(self.runmodel_object.qoi_list[-len(x_run) :]) # Accept the states with g <= g_level - ind_accept = np.where(g_temp <= self.performance_threshold_per_level[step - 1])[0] - for ii in ind_accept: - self.samples[step][(i + 1) * n_keep + ind_false[ii]] = x_run[ii] - self.performance_function_per_level[step][(i + 1) * n_keep + ind_false[ii]] = g_temp[ii] + ind_accept = np.where(response_function_values <= self.performance_threshold_per_level[conditional_level - 1])[0] + for j in ind_accept: + self.samples[conditional_level][(i + 1) * n_keep + ind_false[j]] = x_run[j] + self.performance_function_per_level[conditional_level][(i + 1) * n_keep + ind_false[j]] = response_function_values[j] # Reject the states with g > g_level - ind_reject = np.where(g_temp > self.performance_threshold_per_level[step - 1])[0] - for ii in ind_reject: - self.samples[step][(i + 1) * n_keep + ind_false[ii]] =\ - self.samples[step][i * n_keep + ind_false[ii]] - self.performance_function_per_level[step][(i + 1) * n_keep + ind_false[ii]] = \ - self.performance_function_per_level[step][i * n_keep + ind_false[ii]] + ind_reject = np.where(response_function_values > self.performance_threshold_per_level[conditional_level - 1])[0] + for k in ind_reject: + self.samples[conditional_level][(i + 1) * n_keep + ind_false[k]] =\ + self.samples[conditional_level][i * n_keep + ind_false[k]] + self.performance_function_per_level[conditional_level][(i + 1) * n_keep + ind_false[k]] = \ + self.performance_function_per_level[conditional_level][i * n_keep + ind_false[k]] - g_ind = np.argsort(self.performance_function_per_level[step]) - self.performance_threshold_per_level.append(self.performance_function_per_level[step][g_ind[n_keep]]) + g_ind = np.argsort(self.performance_function_per_level[conditional_level]) + self.performance_threshold_per_level.append(self.performance_function_per_level[conditional_level][g_ind[n_keep]]) # Estimate coefficient of variation of conditional probability of first level - d1, d2 = self._compute_coefficient_of_variation(step) - d12.append(d1 ** 2) - d22.append(d2 ** 2) + independent_intermediate_cov, dependent_intermediate_cov = self._compute_intermediate_cov(conditional_level) + independent_chain_cov_squared.append(independent_intermediate_cov ** 2) + dependent_chain_cov_squared.append(dependent_intermediate_cov ** 2) - self.logger.info("UQpy: Subset Simulation, conditional level " + str(step) + " complete.") + self.logger.info("UQpy: Subset Simulation, conditional level " + str(conditional_level) + " complete.") - n_fail = len([value for value in self.performance_function_per_level[step] if value < 0]) + n_fail = len([value for value in self.performance_function_per_level[conditional_level] if value < 0]) - failure_probability = (self.conditional_probability ** step * n_fail / self.nsamples_per_subset) - probability_cov_independent = np.sqrt(np.sum(d12)) - probability_cov_dependent = np.sqrt(np.sum(d22)) + failure_probability = (self.conditional_probability ** conditional_level * n_fail / self.nsamples_per_subset) + probability_cov_independent = np.sqrt(np.sum(independent_chain_cov_squared)) + probability_cov_dependent = np.sqrt(np.sum(dependent_chain_cov_squared)) return failure_probability, probability_cov_independent, probability_cov_dependent - def _compute_coefficient_of_variation(self, step): - # Here, we assume that the initial samples are drawn to be uncorrelated such that the correction factors do not - # need to be computed. - if step == 0: + def _compute_intermediate_cov(self, conditional_level: int): + """Computes the coefficient of variation of the intermediate failure probability + + Assumes the initial samples are uncorrelated so correction factors do not need to be computed. + + :param conditional_level: Index :math:`i` for the intermediate subset + :return: independent_chains_cov, dependent_chains_cov + """ + if conditional_level == 0: independent_chains_cov = np.sqrt((1 - self.conditional_probability) / (self.conditional_probability * self.nsamples_per_subset)) dependent_chains_cov = np.sqrt((1 - self.conditional_probability) / (self.conditional_probability * self.nsamples_per_subset)) else: - n_c = int(self.conditional_probability * self.nsamples_per_subset) - n_s = int(1 / self.conditional_probability) - indicator = np.reshape(self.performance_function_per_level[step] < self.performance_threshold_per_level[step], (n_s, n_c)) - gamma = self._correlation_factor_gamma(indicator, n_s, n_c) - g_temp = np.reshape(self.performance_function_per_level[step], (n_s, n_c)) - beta_hat = self._correlation_factor_beta(g_temp, step) + n_chains = int(self.conditional_probability * self.nsamples_per_subset) + n_samples_per_chain = int(1 / self.conditional_probability) + indicator = np.reshape(self.performance_function_per_level[conditional_level] < self.performance_threshold_per_level[conditional_level], + (n_samples_per_chain, n_chains)) + gamma = self._correlation_factor_gamma(indicator, n_samples_per_chain, n_chains) + response_function_values = np.reshape(self.performance_function_per_level[conditional_level], (n_samples_per_chain, n_chains)) + beta_hat = self._correlation_factor_beta(response_function_values, conditional_level) independent_chains_cov = np.sqrt(((1 - self.conditional_probability) / (self.conditional_probability * self.nsamples_per_subset)) @@ -229,40 +225,52 @@ def _compute_coefficient_of_variation(self, step): return independent_chains_cov, dependent_chains_cov - # Computes the conventional correlation factor gamma from Au and Beck - def _correlation_factor_gamma(self, indicator, n_s, n_c): - gam = np.zeros(n_s - 1) - r = np.zeros(n_s) + def _correlation_factor_gamma(self, indicator: np.ndarray, n_samples_per_chain: int, n_chains: int): + """Computes the conventional correlation factor :math:`\gamma` as defined by Au and Beck 2001 + + :param indicator: Intermediate indicator function :math:`I_{Conditional Level}(\cdot)` + :param n_samples_per_chain: Number of samples per chain in the MCMC algorithm + :param n_chains: Number of chains in the MCMC algorithm + :return: + """ + gamma = np.zeros(n_samples_per_chain - 1) + r = np.zeros(n_samples_per_chain) ii = indicator * 1 - r_ = ii @ ii.T / n_c - self.conditional_probability ** 2 + r_ = ii @ ii.T / n_chains - self.conditional_probability ** 2 for i in range(r_.shape[0]): r[i] = np.sum(np.diag(r_, i)) / (r_.shape[0] - i) r0 = 0.1 * (1 - 0.1) r = r / r0 - for i in range(n_s - 1): - gam[i] = (1 - ((i + 1) / n_s)) * r[i + 1] - gam = 2 * np.sum(gam) + for i in range(n_samples_per_chain - 1): + gamma[i] = (1 - ((i + 1) / n_samples_per_chain)) * r[i + 1] + gamma = 2 * np.sum(gamma) - return gam + return gamma - # Computes the updated correlation factor beta from Shields et al. - def _correlation_factor_beta(self, g, step): + def _correlation_factor_beta(self, response_function_values, conditional_level: int): + """Computes the updated correlation factor :math:`beta` from Shields, Giovanis, Sundar 2021 + + :param response_function_values: The response function :math:`G` evaluated at the conditional level :math:`i` + :param conditional_level: Index :math:`i` for the intermediate subset + :return: + """ beta = 0 - for i in range(np.shape(g)[1]): - for j in range(i + 1, np.shape(g)[1]): - if g[0, i] == g[0, j]: + for i in range(np.shape(response_function_values)[1]): + for j in range(i + 1, np.shape(response_function_values)[1]): + if response_function_values[0, i] == response_function_values[0, j]: beta += 1 beta *= 2 - ar = np.asarray(self.mcmc_objects[step].acceptance_rate) - ar_mean = np.mean(ar) + acceptance_rate = np.asarray(self.mcmc_objects[conditional_level].acceptance_rate) + mean_acceptance_rate = np.mean(acceptance_rate) - factor = sum((1 - (i + 1) * np.shape(g)[0] / np.shape(g)[1]) * (1 - ar_mean) for i in range(np.shape(g)[0] - 1)) + factor = sum((1 - (i + 1) * np.shape(response_function_values)[0] / np.shape(response_function_values)[1]) * (1 - mean_acceptance_rate) + for i in range(np.shape(response_function_values)[0] - 1)) factor = factor * 2 + 1 - beta = beta / np.shape(g)[1] * factor + beta = beta / np.shape(response_function_values)[1] * factor return beta