From 9cd9d75625223176b5927ba9bca18667a4509bbe Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 12:39:34 -0400 Subject: [PATCH 01/17] initial commit of the InverseFORM class --- .../reliability/taylor_series/InverseFORM.py | 204 ++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 src/UQpy/reliability/taylor_series/InverseFORM.py diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py new file mode 100644 index 00000000..0e282607 --- /dev/null +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -0,0 +1,204 @@ +import logging +import warnings +import numpy as np +import scipy.stats as stats +from typing import Union +from beartype import beartype +from UQpy.distributions import * +from UQpy.transformations import * +from UQpy.run_model.RunModel import RunModel +from UQpy.utilities.ValidationTypes import PositiveInteger +from UQpy.reliability.taylor_series.baseclass.TaylorSeries import TaylorSeries + +warnings.filterwarnings('ignore') + + +class InverseFORM(TaylorSeries): + + @beartype + def __init__( + self, + distributions: Union[None, Distribution, list[Distribution]], + runmodel_object: RunModel, + p_fail: Union[None, float] = 0.05, + beta: Union[None, float] = None, + seed_x: Union[list, np.ndarray] = None, + seed_u: Union[list, np.ndarray] = None, + df_step: Union[int, float] = 0.01, + corr_x: Union[list, np.ndarray] = None, + corr_z: Union[list, np.ndarray] = None, + max_iterations: PositiveInteger = 100, + tolerance_u: Union[float, int, None] = 1e-3, + tolerance_gradient: Union[float, int, None] = 1e-3, + ): + """Class to perform the Inverse First Order Reliability Method. + + Each time :meth:`run` is called the results are appended to the class attributes. + This is a child class of the :class:`TaylorSeries` class. + By default, :code:`tolerance_u` and :code:`tolerance_gradient` must be satisfied for convergence. + Specifying a tolerance as :code:`None` will ignore that criteria, but both cannot be :code:`None`. + + :param distributions: Marginal probability distributions of each random variable. Must be of + type :class:`.DistributionContinuous1D` or :class:`.JointIndependent`. + :param runmodel_object: The computational model. Must be of type :class:`RunModel`. + :param p_fail: Probability of failure defining the feasibility criteria as :math:`||u||=-\\Phi^{-1}(p_{fail})`. + Default: :math:`0.05`. Set to :code:`None` to use :code:`beta` as the feasibility criteria. + Only one of :code:`p_fail` or :code:`beta` may be provided. + :param beta: Hasofer-Lind reliability index defining the feasibility criteria as :math:`||u|| = \\beta`. + Default: :code:`None`. + Only one of code:`p_fail` or :code:`beta` may be provided. + :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` + :param corr_x: Correlation matrix :math:`\mathbf{C_X}` of the random vector :math:`\mathbf{X}` + :param corr_z: Correlation matrix :math:`\mathbf{C_Z}` of the random vector :math:`\mathbf{Z}` + Default: :code:`corr_z` is the identity matrix. + :param max_iterations: Maximum number of iterations for the `HLRF` algorithm. Default: :math:`100` + :param tolerance_u: Convergence threshold for criterion :math:`||\mathbf{U}^{k} - \mathbf{U}^{k-1}||_2 \leq` + :code:`tolerance_u` of the `HLRF` algorithm. + Default: :math:`1.0e-3` + :param tolerance_gradient: Convergence threshold for criterion + :math:`||\\nabla G(\mathbf{U}_i)- \\nabla G(\mathbf{U}_{i-1})||_2 \leq` :code:`tolerance_gradient` + of the `HLRF` algorithm. + Default: :math:`1.0e-3` + """ + self.distributions = distributions + self.runmodel_object = runmodel_object + if (p_fail is not None) and (beta is not None): + raise ValueError('UQpy: Exactly one input (p_fail or beta) must be provided') + elif (p_fail is None) and (beta is None): + raise ValueError('UQpy: Exactly one input (p_fail or beta) must be provided') + elif p_fail is not None: + self.p_fail = p_fail + self.beta = -stats.norm.ppf(self.p_fail) + elif beta is not None: + self.p_fail = stats.norm.cdf(-beta) + self.beta = beta + self.seed_x = seed_x + self.seed_u = seed_u + self.df_step = df_step + self.corr_x = corr_x + self.corr_z = corr_z + self.max_iterations = max_iterations + self.tolerance_u = tolerance_u + self.tolerance_gradient = tolerance_gradient + if (self.tolerance_u is None) and (self.tolerance_gradient is None): + raise ValueError('UQpy: At least one tolerance (tolerance_u or tolerance_gradient) must be provided') + + self.logger = logging.getLogger(__name__) + self.nataf_object = Nataf(distributions=distributions, corr_z=corr_z, corr_x=corr_x) + + # Determine the number of dimensions as the number of random variables + if isinstance(distributions, DistributionContinuous1D): + self.dimension = 1 + elif isinstance(distributions, JointIndependent): + self.dimension = len(distributions.marginals) + elif isinstance(distributions, list): + self.dimension = 0 + for i in range(len(distributions)): + if isinstance(distributions[i], DistributionContinuous1D): + self.dimension += 1 + elif isinstance(distributions[i], JointIndependent): + self.dimension += len(distributions[i].marginals) + + # Initialize attributes + self.alpha: float = np.nan + """Normalized vector in :math:`\\textbf{U}`-space""" + self.alpha_record: list = [] + """Record of alpha :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" + self.beta: float = self.beta + """Feasibility criteria for the optimization :math:`||\\textbf{U}|| = \\beta`""" + self.beta_record: list = [] + """Record of Hasofer-Lind reliability index that defines the feasibility criteria :math:`||u||=\\beta_{HL}`""" + self.design_point_u: list = [] + """Design point in the standard normal space :math:`\\textbf{U}`""" + self.design_point_x: list = [] + """Design point in the parameter space :math:`\\textbf{X}`""" + self.error_record: list = [] + """Record of the final error defined by + :math:`error_u = ||u_new - u||` and :math:`error_{gradient} = || \\nabla G(u_new) - \\nabla G(u)||`""" + self.iteration_record: list = [] + """Record of the number of iterations before algorithm termination""" + self.failure_probability_record: list = [] + """Record of the probability of failure defined by :math:`p_{fail} = \\Phi(-\\beta_{HL})`""" + + def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.ndarray] = None): + """Runs the inverse FORM algorithm. + + :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. + Only one of `seed_x` or `seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. + Only one of `seed_x` or `seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + """ + self.logger.info('UQpy: Running InverseFORM...') + if (seed_x is not None) and (seed_u is not None): + raise ValueError('UQpy: Only one input (seed_x or seed_u) may be provided') + + # Allocate u and the gradient of G(u) as arrays + u = np.zeros([self.max_iterations + 1, self.dimension]) + state_function_gradient = np.zeros([self.max_iterations + 1, self.dimension]) + + # Set initial seed. If both seed_x and seed_u are None, the initial seed is a vector of zeros in U space. + if seed_u is not None: + u[0, :] = seed_u + elif 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 + u[0, :] = Decorrelate(seed_z, self.nataf_object.corr_z).samples_u + + converged = False + iteration = 0 + while (not converged) and (iteration < self.max_iterations): + self.logger.info(f'Number of iteration: {iteration}') + if iteration == 0: + if seed_x is not None: + x = seed_x + else: + seed_z = Correlate(samples_u=u[0, :].reshape(1, -1), corr_z=self.nataf_object.corr_z).samples_z + self.nataf_object.run(samples_z=seed_z.reshape(1, -1), jacobian=True) + x = self.nataf_object.samples_x + else: + z = Correlate(u[iteration, :].reshape(1, -1), self.nataf_object.corr_z).samples_z + self.nataf_object.run(samples_z=z, jacobian=True) + x = self.nataf_object.samples_x + self.logger.info(f'Design Point U: {u[iteration, :]}\nDesign Point X: {x}\n') + state_function_gradient[iteration + 1, :], _, _ = self._derivatives(point_u=u[iteration, :], + point_x=x, + runmodel_object=self.runmodel_object, + nataf_object=self.nataf_object, + df_step=self.df_step, + order='first') + + alpha = state_function_gradient[iteration + 1] + alpha /= np.linalg.norm(state_function_gradient[iteration + 1]) + u[iteration + 1, :] = -alpha * self.beta + + error_u = np.linalg.norm(u[iteration + 1, :] - u[iteration, :]) + error_gradient = np.linalg.norm(state_function_gradient[iteration + 1, :] + - state_function_gradient[iteration, :]) + + converged_u = True if (self.tolerance_u is None) \ + else (error_u <= self.tolerance_u) + converged_gradient = True if (self.tolerance_gradient is None) \ + else (error_gradient <= self.tolerance_gradient) + converged = converged_u and converged_gradient + + if not converged: + iteration += 1 + + if iteration >= self.max_iterations: + self.logger.info(f'UQpy: Maximum number of iterations {self.max_iterations} reached before convergence') + else: + self.alpha_record.append(alpha) + self.beta_record.append(self.beta) + self.design_point_u.append(u[iteration, :]) + self.design_point_x.append(x) + self.error_record.append((error_u, error_gradient)) + self.iteration_record.append(iteration) + self.failure_probability_record.append(self.p_fail) From c3e5459e133c0f338230fee28a6dc08e8eb47d66 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 12:40:26 -0400 Subject: [PATCH 02/17] tests for accuracy and behavior under various inputs for InverseFORM class --- tests/unit_tests/reliability/example_7_2.py | 18 +++ .../reliability/test_inverse_form.py | 106 ++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 tests/unit_tests/reliability/example_7_2.py create mode 100644 tests/unit_tests/reliability/test_inverse_form.py diff --git a/tests/unit_tests/reliability/example_7_2.py b/tests/unit_tests/reliability/example_7_2.py new file mode 100644 index 00000000..cbeb6f33 --- /dev/null +++ b/tests/unit_tests/reliability/example_7_2.py @@ -0,0 +1,18 @@ +import numpy as np + + +def performance_function(samples=None): + """Performance function from Chapter 7 Example 7.2 from Du 2005""" + elastic_modulus = 30e6 + length = 100 + width = 2 + height = 4 + d_0 = 3 + + g = np.zeros(samples.shape[0]) + for i in range(samples.shape[0]): + x = (samples[i, 0] / width**2) ** 2 + y = (samples[i, 1] / height**2) ** 2 + d = ((4 * length**3) / (elastic_modulus * width * height)) * np.sqrt(x + y) + g[i] = d_0 - d + return g diff --git a/tests/unit_tests/reliability/test_inverse_form.py b/tests/unit_tests/reliability/test_inverse_form.py new file mode 100644 index 00000000..3233dde1 --- /dev/null +++ b/tests/unit_tests/reliability/test_inverse_form.py @@ -0,0 +1,106 @@ +import os +import pytest +import numpy as np +from scipy import stats +from UQpy.distributions import Normal +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel +from UQpy.reliability.taylor_series import InverseFORM + + +@pytest.fixture +def inverse_form(): + """Example 7.2 from Chapter 7 of Du 2005 + + Distributions are :math:`P_X \\sim N(500, 100)`, :math:`P_Y \\sim N(1000, 100)` + + Solution from the reference is :math:`u^*=(1.7367, 0.16376)` + """ + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + return InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=0.04054, + tolerance_u=1e-5, + tolerance_gradient=1e-5) + + +def test_no_seed(inverse_form): + """Solution to Chapter 7 Example 7.2 from Du 2005""" + inverse_form.run() + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_seed_x(inverse_form): + seed_x = np.array([625, 900]) + inverse_form.run(seed_x=seed_x) + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_seed_u(inverse_form): + seed_u = np.array([2.4, -1.0]) + inverse_form.run(seed_u=seed_u) + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_both_seeds(inverse_form): + """Expected behavior is to raise ValueError and inform user only one input may be provided""" + seed_x = np.array([1, 2]) + seed_u = np.array([3, 4]) + with pytest.raises(ValueError, match='UQpy: Only one input .* may be provided'): + inverse_form.run(seed_u=seed_u, seed_x=seed_x) + + +def test_neither_tolerance(): + """Expected behavior is to raise ValueError and inform user at least one tolerance must be provided""" + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + with pytest.raises(ValueError, match='UQpy: At least one tolerance .* must be provided'): + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=0.04054, + tolerance_u=None, + tolerance_gradient=None) + + +def test_beta(): + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + # with pytest.raises(ValueError, match='UQpy: At least one tolerance .* must be provided'): + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=None, + beta=-stats.norm.ppf(0.04054)) + inverse_form.run() + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-3) + + +def test_no_beta_no_pfail(): + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + with pytest.raises(ValueError, match='UQpy: Exactly one input .* must be provided'): + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=None, + beta=None) + From 48b6204ba221f441f055c0e68e4b3c5bb70ca3e9 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 12:40:59 -0400 Subject: [PATCH 03/17] documentation for InverseFORM as a child of the TaylorSeries class --- docs/source/reliability/inverse_form.rst | 74 +++++++++++++++++++ docs/source/reliability/taylor_series.rst | 1 + .../reliability/taylor_series/__init__.py | 1 + 3 files changed, 76 insertions(+) create mode 100644 docs/source/reliability/inverse_form.rst diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst new file mode 100644 index 00000000..79bf8af8 --- /dev/null +++ b/docs/source/reliability/inverse_form.rst @@ -0,0 +1,74 @@ +Inverse FORM +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In FORM :cite:`FORM_XDu`, the performance function is linearized according to + +.. math:: G(\textbf{U}) \approx G(\textbf{U}^\star) + \nabla G(\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^\intercal + +where :math:`\textbf{U}^\star` is the expansion point, :math:`G(\textbf{U})` is the performance function evaluated in +the standard normal space and :math:`\nabla G(\textbf{U}^\star)` is the gradient of :math:`G(\textbf{U})` evaluated at +:math:`\textbf{U}^\star`. The probability failure is approximated as + +.. math:: P_{fail} = \Phi(-\beta_{HL}) + +where :math:`\Phi(\cdot)` is the standard normal cumulative distribution function and :math:`\beta_{HL}=||\textbf{U}^*||` +is the norm of the design point known as the Hasofer-Lind reliability index. + +The goal of the inverse FORM algorithm is to find a design point :math:`\textbf{U}^\star` that minimizes the performance +function :math:`G(\textbf{U})` and lies on the hypersphere defined by :math:`||\textbf{U}^*|| = \beta_{HL}`, or +equivalently :math:`||\textbf{U}^*|| = -\Phi^{-1}(p_{fail})`. The default convergence criteria used for this algorithm +are: + +.. math:: \text{tolerance}_{\textbf{U}}:\ ||\textbf{U}_i - \textbf{U}_{i-1}||_2 \leq 10^{-3} +.. math:: \text{tolerance}_{\nabla G(\textbf{U})}:\ ||\nabla G(\textbf{U}_i)- \nabla G(\textbf{U}_{i-1})||_2 \leq 10^{-3} + + +Problem Statement +----- +Compute :math:`u^* = \text{argmin}\ G(\textbf{U})` such that :math:`||\textbf{U}||=\beta`. + +The feasibility criteria :math:`||\textbf{U}||=\beta` may be equivalently defined as +:math:`\beta = -\Phi^{-1}(p_{fail})`, where :math:`\Phi^{-1}(\cdot)` is the inverse standard normal CDF. + +Algorithm +----- +This method implements a gradient descent algorithm to solve the optimization problem within the tolerances specified by +:math:`\text{tolerance}_{\textbf{U}}` (:code:`tolerance_u`) and :math:`\text{tolerance}_{\nabla G(\textbf{U})}` (:code:`tolerance_gradient`). + +0. Define :math:`u_0` and :math:`\beta` directly or :math:`\beta=-\Phi^{-1}(p_\text{fail})` +1. Set :math:`u=u_0` and :math:`\text{converged}=False` +2. While not :math:`\text{converged}`: + a. :math:`\alpha = \frac{\nabla G(u)}{|| \nabla G(u) ||}` + b. :math:`u_{new} = - \alpha \beta` + c. If :math:`||u_{new} - u || \leq \text{tolerance}_u` and/or :math:`|| \nabla G(u_{new}) - \nabla G(u) || \leq \text{tolerance}_{\nabla G(u)}`; + set :math:`\text{converged}=True` + else; + :math:`u = u_{new}` + +The :class:`.InverseFORM` class is imported using the following command: + +>>> from UQpy.reliability.taylor_series import InverseFORM + +Methods +""""""" +.. autoclass:: UQpy.reliability.taylor_series.InverseFORM + :members: run + +Attributes +"""""""""" +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha_record +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta_record +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.design_point_u +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.design_point_x +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.error_record +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.iteration_record +.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.failure_probability_record + + +Examples +"""""""""" +.. toctree:: + + FORM Examples <../auto_examples/reliability/inverse_form/index> \ No newline at end of file diff --git a/docs/source/reliability/taylor_series.rst b/docs/source/reliability/taylor_series.rst index 29ff80e0..9687f597 100644 --- a/docs/source/reliability/taylor_series.rst +++ b/docs/source/reliability/taylor_series.rst @@ -21,6 +21,7 @@ respectively. These classes can be imported in a python script using the followi FORM
SORM + InverseFORM diff --git a/src/UQpy/reliability/taylor_series/__init__.py b/src/UQpy/reliability/taylor_series/__init__.py index c5af5bab..e2797682 100644 --- a/src/UQpy/reliability/taylor_series/__init__.py +++ b/src/UQpy/reliability/taylor_series/__init__.py @@ -1,3 +1,4 @@ from UQpy.reliability.taylor_series.FORM import FORM from UQpy.reliability.taylor_series.SORM import SORM +from UQpy.reliability.taylor_series.InverseFORM import InverseFORM from UQpy.reliability.taylor_series.baseclass.TaylorSeries import TaylorSeries From 55db785a8150a2d84bfb829aa0322b857e99dd0b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 12:41:34 -0400 Subject: [PATCH 04/17] scripts to config and run examples for InverseFORM --- docs/code/reliability/inverse_form/README.rst | 4 + .../inverse_form/inverse_form_cantilever.py | 73 +++++++++++++++++++ .../inverse_form/performance_function.py | 18 +++++ 3 files changed, 95 insertions(+) create mode 100644 docs/code/reliability/inverse_form/README.rst create mode 100644 docs/code/reliability/inverse_form/inverse_form_cantilever.py create mode 100644 docs/code/reliability/inverse_form/performance_function.py diff --git a/docs/code/reliability/inverse_form/README.rst b/docs/code/reliability/inverse_form/README.rst new file mode 100644 index 00000000..b3c07358 --- /dev/null +++ b/docs/code/reliability/inverse_form/README.rst @@ -0,0 +1,4 @@ +Inverse FORM Examples +^^^^^ + +The following example is Example 7.2 from Chapter 7 of :cite:`FORM_XDu`. \ No newline at end of file diff --git a/docs/code/reliability/inverse_form/inverse_form_cantilever.py b/docs/code/reliability/inverse_form/inverse_form_cantilever.py new file mode 100644 index 00000000..ed7b2216 --- /dev/null +++ b/docs/code/reliability/inverse_form/inverse_form_cantilever.py @@ -0,0 +1,73 @@ +""" +Inverse FORM - Cantilever Beam +----- + +A cantilever beam (example 7.2 in :cite:`FORM_XDu`) is considered to fail if the displacement at the tip exceeds the +threshold :math:`D_0`. The performance function :math:`G(\textbf{U})` of this problem is given by + +.. math:: G = D_0 - \frac{4L^3}{Ewt} \sqrt{ \left(\frac{P_x}{w^2}\right)^2 + \left(\frac{P_y}{t^2}\right)^2} + +Where the external forces are modeled as random variables :math:`P_x \sim N(500, 100)` and :math:`P_y \sim N(1000,100)`. +The constants in the problem are length (:math:`L=100`), elastic modulus (:math:E=30\times 10^6), cross section width +(:math:`w=2`) and cross section height (:math:`t=4`). + +""" +# %% md +# +# Import the necessary modules. + +# %% +import numpy as np +import matplotlib.pyplot as plt +plt.style.use('ggplot') +from scipy import stats +from UQpy.distributions import Normal +from UQpy.reliability.taylor_series import InverseFORM +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel + +# %% md +# +# Next, we initialize the :code:`RunModel` object. +# The file defining the performance function file can be found on the UQpy GitHub. +# It contains a function :code:`cantilever_beam` to compute the performance function :math:`G(\textbf{U})`. + +# %% + +model = PythonModel(model_script='performance_function.py', model_object_name="cantilever_beam") +runmodel_object = RunModel(model=model) + +# %% md +# +# Next, we define the external forces in the :math:`x` and :math:`y` direction as distributions that will be passed into +# :code:`FORM`. Along with the distributions, :code:`FORM` takes in the previously defined :code:`runmodel_object`, +# the specified probability of failure, and the tolerances. These tolerances are smaller than the defaults to ensure +# convergence with the level of accuracy given in the problem. + +# %% + +p_fail = 0.04054 +distributions = [Normal(500, 100), Normal(1_000, 100)] +inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=p_fail, + tolerance_u=1e-5, + tolerance_gradient=1e-5) + +# %% md +# +# With everything defined we are ready to run the inverse first-order reliability method and print the results. +# The solution to this problem given by Du is :math:`\textbf{U}^*=(1.7367, 0.16376)` with a reliability index of +# :math:`\beta_{HL}=||\textbf{U}^*||=1.7444` and probability of failure of +# :math:`p_{fail} = \Phi(-\beta_{HL})=\Phi(-1.7444)=0.04054`. We expect this problem to converge in 4 iterations. +# We confirm our design point matches this length, and therefore has a probability of failure specified by our input. + +# %% + +inverse_form.run() +beta = np.linalg.norm(inverse_form.design_point_u) +print('Design point in standard normal space (u^*):', inverse_form.design_point_u[0]) +print('Design point in original space:', inverse_form.design_point_x[0]) +print('Hasofer-Lind reliability index:', beta) +print('Probability of failure at design point:', stats.norm.cdf(-beta)) +print('Number of iterations:', inverse_form.iteration_record[0]) \ No newline at end of file diff --git a/docs/code/reliability/inverse_form/performance_function.py b/docs/code/reliability/inverse_form/performance_function.py new file mode 100644 index 00000000..4bfbcecf --- /dev/null +++ b/docs/code/reliability/inverse_form/performance_function.py @@ -0,0 +1,18 @@ +import numpy as np + + +def cantilever_beam(samples=None): + """Performance function from Chapter 7 Example 7.2 from Du 2005""" + elastic_modulus = 30e6 + length = 100 + width = 2 + height = 4 + d_0 = 3 + + g = np.zeros(samples.shape[0]) + for i in range(samples.shape[0]): + x = (samples[i, 0] / width**2) ** 2 + y = (samples[i, 1] / height**2) ** 2 + d = ((4 * length**3) / (elastic_modulus * width * height)) * np.sqrt(x + y) + g[i] = d_0 - d + return g From bd5c57de1327622eaae61cc6dcceb38e8e08678c Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 14:15:29 -0400 Subject: [PATCH 05/17] added InverseFORM to paths so the examples would show up in the ReadTheDocs --- docs/source/conf.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/source/conf.py b/docs/source/conf.py index b4891bf8..2f4e7635 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -106,6 +106,7 @@ "../code/reliability/form", "../code/reliability/sorm", "../code/reliability/subset_simulation", + "../code/reliability/inverse_form", "../code/surrogates/srom", "../code/surrogates/gpr", "../code/surrogates/pce", @@ -148,6 +149,7 @@ "auto_examples/reliability/form", "auto_examples/reliability/sorm", "auto_examples/reliability/subset_simulation", + "auto_examples/reliability/inverse_form", "auto_examples/surrogates/srom", "auto_examples/surrogates/gpr", "auto_examples/surrogates/pce", From 935d45ddb3cb849aa387e9774561d063832cb85a Mon Sep 17 00:00:00 2001 From: connor-krill Date: Tue, 15 Aug 2023 14:16:14 -0400 Subject: [PATCH 06/17] improved docstrings and markdown for cleaniness and consistency in readthedocs --- docs/code/reliability/inverse_form/README.rst | 4 +--- .../inverse_form/inverse_form_cantilever.py | 22 +++++++++---------- .../{performance_function.py => local_pfn.py} | 6 +++++ docs/source/reliability/inverse_form.rst | 12 +++++----- .../reliability/taylor_series/InverseFORM.py | 7 +++--- 5 files changed, 28 insertions(+), 23 deletions(-) rename docs/code/reliability/inverse_form/{performance_function.py => local_pfn.py} (87%) diff --git a/docs/code/reliability/inverse_form/README.rst b/docs/code/reliability/inverse_form/README.rst index b3c07358..44cdc4d1 100644 --- a/docs/code/reliability/inverse_form/README.rst +++ b/docs/code/reliability/inverse_form/README.rst @@ -1,4 +1,2 @@ Inverse FORM Examples -^^^^^ - -The following example is Example 7.2 from Chapter 7 of :cite:`FORM_XDu`. \ No newline at end of file +^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/code/reliability/inverse_form/inverse_form_cantilever.py b/docs/code/reliability/inverse_form/inverse_form_cantilever.py index ed7b2216..b5bfcb79 100644 --- a/docs/code/reliability/inverse_form/inverse_form_cantilever.py +++ b/docs/code/reliability/inverse_form/inverse_form_cantilever.py @@ -1,25 +1,25 @@ """ Inverse FORM - Cantilever Beam ------ +----------------------------------- +The following example is Example 7.2 from Chapter 7 of :cite:`FORM_XDu`. -A cantilever beam (example 7.2 in :cite:`FORM_XDu`) is considered to fail if the displacement at the tip exceeds the -threshold :math:`D_0`. The performance function :math:`G(\textbf{U})` of this problem is given by +A cantilever beam is considered to fail if the displacement at the tip exceeds the threshold :math:`D_0`. +The performance function :math:`G(\\textbf{U})` of this problem is given by -.. math:: G = D_0 - \frac{4L^3}{Ewt} \sqrt{ \left(\frac{P_x}{w^2}\right)^2 + \left(\frac{P_y}{t^2}\right)^2} +.. math:: G = D_0 - \\frac{4L^3}{Ewt} \\sqrt{ \\left(\\frac{P_x}{w^2}\\right)^2 + \\left(\\frac{P_y}{t^2}\\right)^2} Where the external forces are modeled as random variables :math:`P_x \sim N(500, 100)` and :math:`P_y \sim N(1000,100)`. -The constants in the problem are length (:math:`L=100`), elastic modulus (:math:E=30\times 10^6), cross section width -(:math:`w=2`) and cross section height (:math:`t=4`). +The constants in the problem are length (:math:`L=100`), elastic modulus (:math:`E=30\\times 10^6`), cross section width +(:math:`w=2`), cross section height (:math:`t=4`), and :math:`D_0=3`. """ # %% md # -# Import the necessary modules. +# First, we import the necessary modules. # %% + import numpy as np -import matplotlib.pyplot as plt -plt.style.use('ggplot') from scipy import stats from UQpy.distributions import Normal from UQpy.reliability.taylor_series import InverseFORM @@ -34,7 +34,7 @@ # %% -model = PythonModel(model_script='performance_function.py', model_object_name="cantilever_beam") +model = PythonModel(model_script='local_pfn.py', model_object_name="cantilever_beam") runmodel_object = RunModel(model=model) # %% md @@ -70,4 +70,4 @@ print('Design point in original space:', inverse_form.design_point_x[0]) print('Hasofer-Lind reliability index:', beta) print('Probability of failure at design point:', stats.norm.cdf(-beta)) -print('Number of iterations:', inverse_form.iteration_record[0]) \ No newline at end of file +print('Number of iterations:', inverse_form.iteration_record[0]) diff --git a/docs/code/reliability/inverse_form/performance_function.py b/docs/code/reliability/inverse_form/local_pfn.py similarity index 87% rename from docs/code/reliability/inverse_form/performance_function.py rename to docs/code/reliability/inverse_form/local_pfn.py index 4bfbcecf..4f40a8a7 100644 --- a/docs/code/reliability/inverse_form/performance_function.py +++ b/docs/code/reliability/inverse_form/local_pfn.py @@ -1,3 +1,9 @@ +""" + +Auxiliary file +============================================== + +""" import numpy as np diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst index 79bf8af8..6dde4b65 100644 --- a/docs/source/reliability/inverse_form.rst +++ b/docs/source/reliability/inverse_form.rst @@ -9,7 +9,7 @@ where :math:`\textbf{U}^\star` is the expansion point, :math:`G(\textbf{U})` is the standard normal space and :math:`\nabla G(\textbf{U}^\star)` is the gradient of :math:`G(\textbf{U})` evaluated at :math:`\textbf{U}^\star`. The probability failure is approximated as -.. math:: P_{fail} = \Phi(-\beta_{HL}) +.. math:: p_{fail} = \Phi(-\beta_{HL}) where :math:`\Phi(\cdot)` is the standard normal cumulative distribution function and :math:`\beta_{HL}=||\textbf{U}^*||` is the norm of the design point known as the Hasofer-Lind reliability index. @@ -24,14 +24,14 @@ are: Problem Statement ------ +----------------- Compute :math:`u^* = \text{argmin}\ G(\textbf{U})` such that :math:`||\textbf{U}||=\beta`. The feasibility criteria :math:`||\textbf{U}||=\beta` may be equivalently defined as :math:`\beta = -\Phi^{-1}(p_{fail})`, where :math:`\Phi^{-1}(\cdot)` is the inverse standard normal CDF. Algorithm ------ +----------------- This method implements a gradient descent algorithm to solve the optimization problem within the tolerances specified by :math:`\text{tolerance}_{\textbf{U}}` (:code:`tolerance_u`) and :math:`\text{tolerance}_{\nabla G(\textbf{U})}` (:code:`tolerance_gradient`). @@ -50,12 +50,12 @@ The :class:`.InverseFORM` class is imported using the following command: >>> from UQpy.reliability.taylor_series import InverseFORM Methods -""""""" +----------------- .. autoclass:: UQpy.reliability.taylor_series.InverseFORM :members: run Attributes -"""""""""" +----------------- .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha_record .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta @@ -71,4 +71,4 @@ Examples """""""""" .. toctree:: - FORM Examples <../auto_examples/reliability/inverse_form/index> \ No newline at end of file + InverseFORM Examples <../auto_examples/reliability/inverse_form/index> \ No newline at end of file diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index 0e282607..6ec2adff 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -111,16 +111,17 @@ def __init__( self.alpha_record: list = [] """Record of alpha :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" self.beta: float = self.beta - """Feasibility criteria for the optimization :math:`||\\textbf{U}|| = \\beta`""" + """Feasibility criteria for the optimization :math:`||\\textbf{U}|| = \\beta_{HL}`""" self.beta_record: list = [] - """Record of Hasofer-Lind reliability index that defines the feasibility criteria :math:`||u||=\\beta_{HL}`""" + """Record of Hasofer-Lind reliability index that defines the feasibility criteria + :math:`||\\textbf{U}||=\\beta_{HL}`""" self.design_point_u: list = [] """Design point in the standard normal space :math:`\\textbf{U}`""" self.design_point_x: list = [] """Design point in the parameter space :math:`\\textbf{X}`""" self.error_record: list = [] """Record of the final error defined by - :math:`error_u = ||u_new - u||` and :math:`error_{gradient} = || \\nabla G(u_new) - \\nabla G(u)||`""" + :math:`error_u = ||u_{new} - u||` and :math:`error_{\\nabla G(u)} = || \\nabla G(u_{new}) - \\nabla G(u)||`""" self.iteration_record: list = [] """Record of the number of iterations before algorithm termination""" self.failure_probability_record: list = [] From b74b418e173572f403ddb13e4b6d1c32eb4204d5 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 16 Aug 2023 10:06:22 -0400 Subject: [PATCH 07/17] corrected vector notation in math statement --- docs/source/reliability/index.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/reliability/index.rst b/docs/source/reliability/index.rst index d3010039..0464cbed 100644 --- a/docs/source/reliability/index.rst +++ b/docs/source/reliability/index.rst @@ -3,7 +3,7 @@ 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 +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}` 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 @@ -22,8 +22,8 @@ transformation and can support reliability analysis for problems with arbitraril This module contains functionality for all reliability methods supported in :py:mod:`UQpy`. The module currently contains the following classes: -- :class:`.TaylorSeries`: Class to perform reliability analysis using First Order reliability Method (:class:`FORM`) and Second Order - Reliability Method (:class:`SORM`). +- :class:`.TaylorSeries`: Class to perform reliability analysis using First Order Reliability Method (:class:`FORM`), + Inverse First Order Reliability Method (:class:`InverseFORM`) and Second Order Reliability Method (:class:`SORM`). - :class:`.SubsetSimulation`: Class to perform reliability analysis using subset simulation. From 770256bd9e49f6b65cabcee7747934868c60fd62 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 16 Aug 2023 10:07:09 -0400 Subject: [PATCH 08/17] edited markdown for clarity and notational consistency --- .../inverse_form/inverse_form_cantilever.py | 2 +- docs/source/reliability/inverse_form.rst | 13 ++++++++----- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/docs/code/reliability/inverse_form/inverse_form_cantilever.py b/docs/code/reliability/inverse_form/inverse_form_cantilever.py index b5bfcb79..54f355e7 100644 --- a/docs/code/reliability/inverse_form/inverse_form_cantilever.py +++ b/docs/code/reliability/inverse_form/inverse_form_cantilever.py @@ -59,7 +59,7 @@ # With everything defined we are ready to run the inverse first-order reliability method and print the results. # The solution to this problem given by Du is :math:`\textbf{U}^*=(1.7367, 0.16376)` with a reliability index of # :math:`\beta_{HL}=||\textbf{U}^*||=1.7444` and probability of failure of -# :math:`p_{fail} = \Phi(-\beta_{HL})=\Phi(-1.7444)=0.04054`. We expect this problem to converge in 4 iterations. +# :math:`p_{fail} = \Phi(-\beta_{HL})=0.04054`. We expect this problem to converge in 4 iterations. # We confirm our design point matches this length, and therefore has a probability of failure specified by our input. # %% diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst index 6dde4b65..f3470823 100644 --- a/docs/source/reliability/inverse_form.rst +++ b/docs/source/reliability/inverse_form.rst @@ -31,9 +31,9 @@ The feasibility criteria :math:`||\textbf{U}||=\beta` may be equivalently define :math:`\beta = -\Phi^{-1}(p_{fail})`, where :math:`\Phi^{-1}(\cdot)` is the inverse standard normal CDF. Algorithm ------------------ +--------- This method implements a gradient descent algorithm to solve the optimization problem within the tolerances specified by -:math:`\text{tolerance}_{\textbf{U}}` (:code:`tolerance_u`) and :math:`\text{tolerance}_{\nabla G(\textbf{U})}` (:code:`tolerance_gradient`). +:code:`tolerance_u` and :code:`tolerance_gradient`. 0. Define :math:`u_0` and :math:`\beta` directly or :math:`\beta=-\Phi^{-1}(p_\text{fail})` 1. Set :math:`u=u_0` and :math:`\text{converged}=False` @@ -50,12 +50,14 @@ The :class:`.InverseFORM` class is imported using the following command: >>> from UQpy.reliability.taylor_series import InverseFORM Methods ------------------ +------- + .. autoclass:: UQpy.reliability.taylor_series.InverseFORM :members: run Attributes ------------------ +---------- + .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha_record .. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta @@ -68,7 +70,8 @@ Attributes Examples -"""""""""" +-------- + .. toctree:: InverseFORM Examples <../auto_examples/reliability/inverse_form/index> \ No newline at end of file From fa0230243a164d895c95a3a051fc3113cda78b4a Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 16 Aug 2023 10:07:46 -0400 Subject: [PATCH 09/17] editted docstrings for clarity and notational consistency --- .../reliability/taylor_series/InverseFORM.py | 29 ++++++++++--------- .../reliability/test_inverse_form.py | 9 +++--- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index 6ec2adff..be4610cf 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -34,9 +34,9 @@ def __init__( """Class to perform the Inverse First Order Reliability Method. Each time :meth:`run` is called the results are appended to the class attributes. - This is a child class of the :class:`TaylorSeries` class. - By default, :code:`tolerance_u` and :code:`tolerance_gradient` must be satisfied for convergence. - Specifying a tolerance as :code:`None` will ignore that criteria, but both cannot be :code:`None`. + By default, :code:`tolerance_u` and :code:`tolerance_gradient` must be satisfied for convergence. + Specifying a tolerance as :code:`None` will ignore that criteria, but both cannot be :code:`None`. + This is a child class of the :class:`TaylorSeries` class. :param distributions: Marginal probability distributions of each random variable. Must be of type :class:`.DistributionContinuous1D` or :class:`.JointIndependent`. @@ -46,19 +46,19 @@ def __init__( Only one of :code:`p_fail` or :code:`beta` may be provided. :param beta: Hasofer-Lind reliability index defining the feasibility criteria as :math:`||u|| = \\beta`. Default: :code:`None`. - Only one of code:`p_fail` or :code:`beta` may be provided. + Only one of :code:`p_fail` or :code:`beta` may be provided. :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. Only one of :code:`seed_x` or :code:`seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. Only one of :code:`seed_x` or :code:`seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` :param corr_x: Correlation matrix :math:`\mathbf{C_X}` of the random vector :math:`\mathbf{X}` :param corr_z: Correlation matrix :math:`\mathbf{C_Z}` of the random vector :math:`\mathbf{Z}` Default: :code:`corr_z` is the identity matrix. :param max_iterations: Maximum number of iterations for the `HLRF` algorithm. Default: :math:`100` - :param tolerance_u: Convergence threshold for criterion :math:`||\mathbf{U}^{k} - \mathbf{U}^{k-1}||_2 \leq` + :param tolerance_u: Convergence threshold for criterion :math:`||\mathbf{U}_i - \mathbf{U}_{i-1}||_2 \leq` :code:`tolerance_u` of the `HLRF` algorithm. Default: :math:`1.0e-3` :param tolerance_gradient: Convergence threshold for criterion @@ -107,11 +107,12 @@ def __init__( # Initialize attributes self.alpha: float = np.nan - """Normalized vector in :math:`\\textbf{U}`-space""" + """Normalized gradient vector in :math:`\\textbf{U}` space""" self.alpha_record: list = [] - """Record of alpha :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" + """Record of :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" self.beta: float = self.beta - """Feasibility criteria for the optimization :math:`||\\textbf{U}|| = \\beta_{HL}`""" + """Hasofer-Lind reliability index that defines the feasibility criteria as + :math:`||\\textbf{U}|| = \\beta_{HL}`""" self.beta_record: list = [] """Record of Hasofer-Lind reliability index that defines the feasibility criteria :math:`||\\textbf{U}||=\\beta_{HL}`""" @@ -131,11 +132,11 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda """Runs the inverse FORM algorithm. :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. - Only one of `seed_x` or `seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. - Only one of `seed_x` or `seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is used as the seed. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. """ self.logger.info('UQpy: Running InverseFORM...') if (seed_x is not None) and (seed_u is not None): diff --git a/tests/unit_tests/reliability/test_inverse_form.py b/tests/unit_tests/reliability/test_inverse_form.py index 3233dde1..1c0e753e 100644 --- a/tests/unit_tests/reliability/test_inverse_form.py +++ b/tests/unit_tests/reliability/test_inverse_form.py @@ -10,11 +10,11 @@ @pytest.fixture def inverse_form(): - """Example 7.2 from Chapter 7 of Du 2005 + """Example 7.2 from Chapter 7 of X. Du 2005 Distributions are :math:`P_X \\sim N(500, 100)`, :math:`P_Y \\sim N(1000, 100)` - - Solution from the reference is :math:`u^*=(1.7367, 0.16376)` + Solution from the reference is :math:`u^*=(1.7367, 0.16376)`. + Tolerances of :math:`1e-5` are used to ensure convergence to level of precision given by Du. """ path = os.path.abspath(os.path.dirname(__file__)) os.chdir(path) @@ -31,7 +31,6 @@ def inverse_form(): def test_no_seed(inverse_form): - """Solution to Chapter 7 Example 7.2 from Du 2005""" inverse_form.run() assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) @@ -81,7 +80,6 @@ def test_beta(): delete_files=True) runmodel_object = RunModel(model=python_model) distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] - # with pytest.raises(ValueError, match='UQpy: At least one tolerance .* must be provided'): inverse_form = InverseFORM(distributions=distributions, runmodel_object=runmodel_object, p_fail=None, @@ -91,6 +89,7 @@ def test_beta(): def test_no_beta_no_pfail(): + """Expected behavior is to raise ValueError and inform the user exactly one in put must be provided""" path = os.path.abspath(os.path.dirname(__file__)) os.chdir(path) python_model = PythonModel(model_script='example_7_2.py', From adafcdfb17e86098adfb8de5db06631cdcf6500b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 16 Aug 2023 10:08:20 -0400 Subject: [PATCH 10/17] added references to the new child class InverseFORM --- docs/source/reliability/taylor_series.rst | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/docs/source/reliability/taylor_series.rst b/docs/source/reliability/taylor_series.rst index 9687f597..78ec77a6 100644 --- a/docs/source/reliability/taylor_series.rst +++ b/docs/source/reliability/taylor_series.rst @@ -1,27 +1,30 @@ Taylor Series ------------- -:class:`.TaylorSeries` is a class that calculates the reliability of a model using the First Order Reliability Method (FORM) -or the Second Order Reliability Method (SORM) based on the first-order and second-order Taylor series expansion -approximation of the performance function, respectively (:cite:`TaylorSeries1`, :cite:`TaylorSeries2`). +:class:`.TaylorSeries` is a class that calculates the reliability of a model using the First Order Reliability Method +(FORM), Inverse First Order Reliability Method (InverseFORM) or Second Order Reliability Method (SORM) based on the +first-order and second-order Taylor series expansion approximation of the performance function, respectively +(:cite:`TaylorSeries1`, :cite:`FORM_XDu`, :cite:`TaylorSeries2`). .. image:: ../_static/Reliability_FORM.png :scale: 40 % :alt: Graphical representation of the FORM. :align: center -The :class:`.TaylorSeries` class is the parent class of the :class:`.FORM` and :class:`.SORM` classes that perform the FORM and SORM, -respectively. These classes can be imported in a python script using the following command: +The :class:`.TaylorSeries` class is the parent class of the :class:`.FORM`, :class:`InverseFORM`, and :class:`.SORM` +classes that perform the FORM, InverseFORM, and SORM, respectively. +These classes can be imported in a python script using the following command: ->>> from UQpy.reliability.taylor_series import FORM, SORM +>>> from UQpy.reliability.taylor_series import FORM, InverseFORM, SORM .. toctree:: :maxdepth: 1 FORM - SORM InverseFORM + SORM + From 5a6a1ce2f61c400bef68bb4c77f01b8633381725 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Thu, 17 Aug 2023 13:30:14 -0400 Subject: [PATCH 11/17] added state_function attribute, correct attribute behavior when max number of iterations is reached --- .../reliability/taylor_series/InverseFORM.py | 33 +++++++++++-------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index be4610cf..991d5c10 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -127,6 +127,8 @@ def __init__( """Record of the number of iterations before algorithm termination""" self.failure_probability_record: list = [] """Record of the probability of failure defined by :math:`p_{fail} = \\Phi(-\\beta_{HL})`""" + self.state_function: list = [] + """State function :math:`G(u)` evaluated at each step in the optimization""" def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.ndarray] = None): """Runs the inverse FORM algorithm. @@ -144,6 +146,7 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda # Allocate u and the gradient of G(u) as arrays u = np.zeros([self.max_iterations + 1, self.dimension]) + state_function = np.full(self.max_iterations + 1, np.nan) state_function_gradient = np.zeros([self.max_iterations + 1, self.dimension]) # Set initial seed. If both seed_x and seed_u are None, the initial seed is a vector of zeros in U space. @@ -170,12 +173,14 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda self.nataf_object.run(samples_z=z, jacobian=True) x = self.nataf_object.samples_x self.logger.info(f'Design Point U: {u[iteration, :]}\nDesign Point X: {x}\n') - state_function_gradient[iteration + 1, :], _, _ = self._derivatives(point_u=u[iteration, :], - point_x=x, - runmodel_object=self.runmodel_object, - nataf_object=self.nataf_object, - df_step=self.df_step, - order='first') + state_function_gradient[iteration + 1, :], qoi, _ = self._derivatives(point_u=u[iteration, :], + point_x=x, + runmodel_object=self.runmodel_object, + nataf_object=self.nataf_object, + df_step=self.df_step, + order='first') + self.logger.info(f'State Function: {qoi}') + self.state_function[iteration] = qoi alpha = state_function_gradient[iteration + 1] alpha /= np.linalg.norm(state_function_gradient[iteration + 1]) @@ -196,11 +201,11 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda if iteration >= self.max_iterations: self.logger.info(f'UQpy: Maximum number of iterations {self.max_iterations} reached before convergence') - else: - self.alpha_record.append(alpha) - self.beta_record.append(self.beta) - self.design_point_u.append(u[iteration, :]) - self.design_point_x.append(x) - self.error_record.append((error_u, error_gradient)) - self.iteration_record.append(iteration) - self.failure_probability_record.append(self.p_fail) + self.alpha_record.append(alpha) + self.beta_record.append(self.beta) + self.design_point_u.append(u[iteration, :]) + self.design_point_x.append(x) + self.error_record.append((error_u, error_gradient)) + self.iteration_record.append(iteration) + self.failure_probability_record.append(self.p_fail) + self.state_function.append(state_function) From d9923886b01190ee6b20b6e280fdba63548ec575 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Thu, 17 Aug 2023 13:34:49 -0400 Subject: [PATCH 12/17] fixed minor bug with state_function attribute --- src/UQpy/reliability/taylor_series/InverseFORM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index 991d5c10..bac070f3 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -180,7 +180,7 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda df_step=self.df_step, order='first') self.logger.info(f'State Function: {qoi}') - self.state_function[iteration] = qoi + state_function[iteration + 1] = qoi alpha = state_function_gradient[iteration + 1] alpha /= np.linalg.norm(state_function_gradient[iteration + 1]) From 738e25b437037f04200b980465b6f3d03e1d5660 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 18 Aug 2023 16:03:27 -0400 Subject: [PATCH 13/17] changed title font to match existing documentation --- docs/source/reliability/inverse_form.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst index f3470823..502b3c9c 100644 --- a/docs/source/reliability/inverse_form.rst +++ b/docs/source/reliability/inverse_form.rst @@ -23,15 +23,15 @@ are: .. math:: \text{tolerance}_{\nabla G(\textbf{U})}:\ ||\nabla G(\textbf{U}_i)- \nabla G(\textbf{U}_{i-1})||_2 \leq 10^{-3} -Problem Statement ------------------ +**Problem Statement** + Compute :math:`u^* = \text{argmin}\ G(\textbf{U})` such that :math:`||\textbf{U}||=\beta`. The feasibility criteria :math:`||\textbf{U}||=\beta` may be equivalently defined as :math:`\beta = -\Phi^{-1}(p_{fail})`, where :math:`\Phi^{-1}(\cdot)` is the inverse standard normal CDF. -Algorithm ---------- +**Algorithm** + This method implements a gradient descent algorithm to solve the optimization problem within the tolerances specified by :code:`tolerance_u` and :code:`tolerance_gradient`. From 6236297eb4c12fc4e38e357b5186b64d20d4f484 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 1 Sep 2023 12:53:06 -0400 Subject: [PATCH 14/17] added logic to prevent seed_x and see_u from being specified at the same time in __init__ --- src/UQpy/reliability/taylor_series/FORM.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index c10c895c..d64300fe 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -133,6 +133,8 @@ def __init__( self.x_record: list = [] """Record of all iteration points in the parameter space **X**.""" + if (seed_x is not None) and (seed_u is not None): + raise ValueError('UQpy: Only one input (seed_x or seed_u) may be provided') if self.seed_u is not None: self.run(seed_u=self.seed_u) elif self.seed_x is not None: From 9704aa0294266d7aaf1b144e697e4be6f195a318 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 1 Sep 2023 12:53:29 -0400 Subject: [PATCH 15/17] added logic to run InverseFORM.run() if seed is provided in init --- src/UQpy/reliability/taylor_series/InverseFORM.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index bac070f3..288053e8 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -49,10 +49,10 @@ def __init__( Only one of :code:`p_fail` or :code:`beta` may be provided. :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. Only one of :code:`seed_x` or :code:`seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. + If either :code:`seed_u` or :code:`seed_x` is provided, then the :py:meth:`run` method will be executed automatically. :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. Only one of :code:`seed_x` or :code:`seed_u` may be provided. - If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. + If either :code:`seed_u` or :code:`seed_x` is provided, then the :py:meth:`run` method will be executed automatically. :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` :param corr_x: Correlation matrix :math:`\mathbf{C_X}` of the random vector :math:`\mathbf{X}` :param corr_z: Correlation matrix :math:`\mathbf{C_Z}` of the random vector :math:`\mathbf{Z}` @@ -130,6 +130,13 @@ def __init__( self.state_function: list = [] """State function :math:`G(u)` evaluated at each step in the optimization""" + if (seed_x is not None) and (seed_u is not None): + raise ValueError('UQpy: Only one input (seed_x or seed_u) may be provided') + if self.seed_u is not None: + self.run(seed_u=self.seed_u) + elif self.seed_x is not None: + self.run(seed_x=self.seed_x) + def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.ndarray] = None): """Runs the inverse FORM algorithm. From 1aa56e07880f308315646d96e262d9aec2a1e29e Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 1 Sep 2023 13:04:25 -0400 Subject: [PATCH 16/17] fixes dummy syntax --- src/UQpy/reliability/taylor_series/InverseFORM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index 288053e8..5ca737d7 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -76,7 +76,7 @@ def __init__( self.p_fail = p_fail self.beta = -stats.norm.ppf(self.p_fail) elif beta is not None: - self.p_fail = stats.norm.cdf(-beta) + self.p_fail = stats.norm.cdf(-1*beta) self.beta = beta self.seed_x = seed_x self.seed_u = seed_u From d8c3b73b40f6334bedb813f466b9938b2fdc5043 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 1 Sep 2023 13:07:01 -0400 Subject: [PATCH 17/17] deleted unneeded attribute --- src/UQpy/reliability/taylor_series/InverseFORM.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py index 5ca737d7..1e4f4fa9 100644 --- a/src/UQpy/reliability/taylor_series/InverseFORM.py +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -110,9 +110,6 @@ def __init__( """Normalized gradient vector in :math:`\\textbf{U}` space""" self.alpha_record: list = [] """Record of :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" - self.beta: float = self.beta - """Hasofer-Lind reliability index that defines the feasibility criteria as - :math:`||\\textbf{U}|| = \\beta_{HL}`""" self.beta_record: list = [] """Record of Hasofer-Lind reliability index that defines the feasibility criteria :math:`||\\textbf{U}||=\\beta_{HL}`"""