diff --git a/docs/code/reliability/inverse_form/README.rst b/docs/code/reliability/inverse_form/README.rst
new file mode 100644
index 000000000..44cdc4d14
--- /dev/null
+++ b/docs/code/reliability/inverse_form/README.rst
@@ -0,0 +1,2 @@
+Inverse FORM Examples
+^^^^^^^^^^^^^^^^^^^^^^^^^
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 000000000..54f355e72
--- /dev/null
+++ b/docs/code/reliability/inverse_form/inverse_form_cantilever.py
@@ -0,0 +1,73 @@
+"""
+Inverse FORM - Cantilever Beam
+-----------------------------------
+The following example is Example 7.2 from Chapter 7 of :cite:`FORM_XDu`.
+
+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}
+
+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`), cross section height (:math:`t=4`), and :math:`D_0=3`.
+
+"""
+# %% md
+#
+# First, we import the necessary modules.
+
+# %%
+
+import numpy as np
+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='local_pfn.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})=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])
diff --git a/docs/code/reliability/inverse_form/local_pfn.py b/docs/code/reliability/inverse_form/local_pfn.py
new file mode 100644
index 000000000..4f40a8a79
--- /dev/null
+++ b/docs/code/reliability/inverse_form/local_pfn.py
@@ -0,0 +1,24 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+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
diff --git a/docs/source/conf.py b/docs/source/conf.py
index b4891bf87..2f4e76353 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",
diff --git a/docs/source/reliability/index.rst b/docs/source/reliability/index.rst
index 51362cfb1..1a6666cd7 100644
--- a/docs/source/reliability/index.rst
+++ b/docs/source/reliability/index.rst
@@ -2,10 +2,10 @@ Reliability
===========
Reliability of a system refers to the assessment of its probability of failure (i.e the system no longer satisfies some
-performance measure), given the model uncertainty in the structural, environmental and load parameters. Given a vector
+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]^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:`\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
.. 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}
@@ -23,8 +23,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.
diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst
new file mode 100644
index 000000000..502b3c9c5
--- /dev/null
+++ b/docs/source/reliability/inverse_form.rst
@@ -0,0 +1,77 @@
+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
+: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`
+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::
+
+ InverseFORM 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 29ff80e03..78ec77a60 100644
--- a/docs/source/reliability/taylor_series.rst
+++ b/docs/source/reliability/taylor_series.rst
@@ -1,27 +1,31 @@
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 <form>
+ InverseFORM <inverse_form>
SORM <sorm>
+
diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py
index acd94e10b..8c167e526 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:
diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py
new file mode 100644
index 000000000..1e4f4fa95
--- /dev/null
+++ b/src/UQpy/reliability/taylor_series/InverseFORM.py
@@ -0,0 +1,215 @@
+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.
+ 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`.
+ :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 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 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}`
+ 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}_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
+ :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(-1*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 gradient vector in :math:`\\textbf{U}` space"""
+ self.alpha_record: list = []
+ """Record of :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`"""
+ self.beta_record: list = []
+ """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_{\\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 = []
+ """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"""
+
+ 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.
+
+ :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.
+ :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.
+ """
+ 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 = 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.
+ 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, :], 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}')
+ state_function[iteration + 1] = qoi
+
+ 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')
+ 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)
diff --git a/src/UQpy/reliability/taylor_series/__init__.py b/src/UQpy/reliability/taylor_series/__init__.py
index c5af5bab2..e27976823 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
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 000000000..cbeb6f335
--- /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 000000000..1c0e753e2
--- /dev/null
+++ b/tests/unit_tests/reliability/test_inverse_form.py
@@ -0,0 +1,105 @@
+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 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)`.
+ 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)
+ 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):
+ 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)]
+ 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():
+ """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',
+ 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)
+