Merge pull request #227 from SURGroup/Development

Development

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  |
 +-----------------------+------------------------------------------------------------------+

docs/code/reliability/form/plot_FORM_linear_function_3d.py → 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)
 

docs/code/reliability/inverse_form/README.rst

@@ -0,0 +1,2 @@
+Inverse FORM Examples
+^^^^^^^^^^^^^^^^^^^^^^^^^

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

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

docs/code/reliability/sorm/plot_SORM_nonlinear_function.py → docs/code/reliability/sorm/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)
 

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

docs/code/reliability/sorm/local_pfn.py

@@ -6,33 +6,31 @@
 """
 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
     d = 2
     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
+
+
+

docs/code/surrogates/pce/plot_pce_sparsity_lars.py

@@ -15,16 +15,16 @@
 # %%
 
 import numpy as np
-import math
-import numpy as np
+
 from UQpy.distributions import Uniform, JointIndependent
 from UQpy.surrogates import *
 
 
 # %% 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`

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",

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  |
 +-----------------------+------------------------------------------------------------------+

docs/source/reliability/index.rst

@@ -3,27 +3,28 @@ 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
 
 
-.. 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.
 
 
 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.