From 00b91193d1d3dfda355a830da5e3071e9a435c99 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 11:37:24 -0400 Subject: [PATCH 01/22] optimized import statements --- .../grassmann/plot_grassmann_distances.py | 9 ++------- .../grassmann/plot_grassmann_karcher.py | 8 +++----- .../grassmann/plot_grassmann_kernel.py | 9 +++------ .../grassmann/plot_grassmann_log_exp.py | 8 +++----- 4 files changed, 11 insertions(+), 23 deletions(-) diff --git a/docs/code/dimension_reduction/grassmann/plot_grassmann_distances.py b/docs/code/dimension_reduction/grassmann/plot_grassmann_distances.py index 9ebe4afb3..05ce844de 100644 --- a/docs/code/dimension_reduction/grassmann/plot_grassmann_distances.py +++ b/docs/code/dimension_reduction/grassmann/plot_grassmann_distances.py @@ -13,19 +13,14 @@ # %% -import numpy as np -import matplotlib import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable -from UQpy import SVDProjection -import sys +import numpy as np +from UQpy import SVDProjection from UQpy.utilities import GrassmannPoint from UQpy.utilities.distances.baseclass.GrassmannianDistance import GrassmannianDistance from UQpy.utilities.distances.grassmannian_distances.GeodesicDistance import GeodesicDistance -from UQpy.dimension_reduction import GrassmannOperations - # %% md # # Generate four random matrices with reduced rank corresponding to the different samples. The samples are stored in diff --git a/docs/code/dimension_reduction/grassmann/plot_grassmann_karcher.py b/docs/code/dimension_reduction/grassmann/plot_grassmann_karcher.py index 05b4de660..3c8ebd3c5 100644 --- a/docs/code/dimension_reduction/grassmann/plot_grassmann_karcher.py +++ b/docs/code/dimension_reduction/grassmann/plot_grassmann_karcher.py @@ -13,13 +13,11 @@ #%% -import numpy as np -import matplotlib import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable -import sys -from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection +import numpy as np + from UQpy.dimension_reduction import GrassmannOperations +from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection #%% md # diff --git a/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py b/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py index baea30bac..30b8d6f57 100644 --- a/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py +++ b/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py @@ -13,15 +13,12 @@ # %% -import numpy as np -import matplotlib import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable +import numpy as np + from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection -from UQpy.dimension_reduction import GrassmannOperations from UQpy.utilities import GrassmannPoint -from UQpy.utilities.kernels import Kernel, ProjectionKernel -import sys +from UQpy.utilities.kernels import ProjectionKernel # %% md # diff --git a/docs/code/dimension_reduction/grassmann/plot_grassmann_log_exp.py b/docs/code/dimension_reduction/grassmann/plot_grassmann_log_exp.py index 5440f68cc..a8ae8927c 100644 --- a/docs/code/dimension_reduction/grassmann/plot_grassmann_log_exp.py +++ b/docs/code/dimension_reduction/grassmann/plot_grassmann_log_exp.py @@ -13,13 +13,11 @@ #%% -import numpy as np -import matplotlib import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable -from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection +import numpy as np + from UQpy.dimension_reduction import GrassmannOperations -import sys +from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection #%% md # From d41643475cbef2e6a73f9c783ce553baf3e07591 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 11:50:04 -0400 Subject: [PATCH 02/22] removed incorrected and unneeded example statements --- .../grassmann/grassmann_operations.rst | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/docs/source/dimension_reduction/grassmann/grassmann_operations.rst b/docs/source/dimension_reduction/grassmann/grassmann_operations.rst index 439b8f2be..7a9630bdf 100644 --- a/docs/source/dimension_reduction/grassmann/grassmann_operations.rst +++ b/docs/source/dimension_reduction/grassmann/grassmann_operations.rst @@ -46,10 +46,6 @@ To use the :meth:`.exp_map` method, one needs to import the :class:`.GrassmannOp >>> from UQpy.dimension_reduction.grassmann_manifold import GrassmannOperations -It can then be used as follows: - ->>> GrassmannOperations.exp_map() - Since :meth:`.exp_map` is a static method, it does not require instantiation of the :class:`.GrassmannOperations` class. .. automethod:: UQpy.dimension_reduction.grassmann_manifold.GrassmannOperations.exp_map @@ -74,10 +70,6 @@ To use the :meth:`.log_map` method, one needs to import the :class:`.GrassmannOp >>> from UQpy.dimension_reduction.grassmann_manifold import GrassmannOperations -It can then be used as follows: - ->>> GrassmannOperations.log_map() - Since :meth:`.log_map` is a static method, it does not require instantiation of the :class:`.GrassmannOperations` class. @@ -98,10 +90,6 @@ To use the :meth:`.karcher_mean` method, one needs to import the :class:`.Grassm >>> from UQpy.dimension_reduction.grassmann_manifold import GrassmannOperations -It can then be used as follows: - ->>> Grassmann.karcher_mean() - Since :meth:`.karcher_mean` is a static method, it does not require instantiation of the :class:`.GrassmannOperations` class. @@ -126,10 +114,6 @@ To use the :meth:`.frechet_variance` method, one needs to import the :class:`.Gr >>> from UQpy.dimension_reduction.grassmann_manifold import GrassmannOperations -It can then be used as follows: - ->>> GrassmannOperations.frechet_variance() - Since :meth:`.frechet_variance` is a static method, it does not require instantiation of the :class:`.GrassmannOperations` class. From 655d0f605456e9cc2379c884fd1f745ecfe8038b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 11:50:30 -0400 Subject: [PATCH 03/22] improved docstrings --- .../grassmann_manifold/GrassmannInterpolation.py | 4 ++-- .../grassmann_manifold/GrassmannOperations.py | 12 ++++++++---- src/UQpy/reliability/taylor_series/FORM.py | 1 + 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannInterpolation.py b/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannInterpolation.py index 44c5b934b..630922ebf 100644 --- a/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannInterpolation.py +++ b/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannInterpolation.py @@ -24,10 +24,10 @@ def __init__(self, interpolation_method: Union[Surrogate, callable, None], :param interpolation_method: Type of interpolation to perform. This may be specified as a :class:`Surrogate` object or a callable function. If :any:`None`, then multi-linear interpolation is performed. :param manifold_data: Data points on the Grassmann manifold. - :param optimization_method: Optimization method for calculating the Karcher mean. See - :py:meth:`.GrassmannOperations.karcher_mean`. :param coordinates: Nodes of the interpolant. :param distance: Distance measure. + :param optimization_method: Optimization method for calculating the Karcher mean. See + :py:meth:`.GrassmannOperations.karcher_mean`. """ self.interpolation_method = interpolation_method diff --git a/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannOperations.py b/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannOperations.py index 78f100aed..f8dfdbe97 100644 --- a/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannOperations.py +++ b/src/UQpy/dimension_reduction/grassmann_manifold/GrassmannOperations.py @@ -15,7 +15,8 @@ class GrassmannOperations: @beartype def __init__(self, grassmann_points: Union[list[Numpy2DFloatArrayOrthonormal], list[GrassmannPoint]], kernel: GrassmannianKernel = ProjectionKernel(), - p: Union[int, str] = "max", optimization_method: str = "GradientDescent", + p: Union[int, str] = "max", + optimization_method: str = "GradientDescent", distance: GrassmannianDistance = GeodesicDistance()): """ The :class:`.GrassmannOperations` class can used in two ways. In the first case, the user can invoke the @@ -49,7 +50,8 @@ def calculate_kernel_matrix(grassmann_points: Union[list[Numpy2DFloatArrayOrthon @beartype def log_map(grassmann_points: Union[list[Numpy2DFloatArrayOrthonormal], list[GrassmannPoint]], reference_point: Union[Numpy2DFloatArrayOrthonormal, GrassmannPoint]) -> list[Numpy2DFloatArray]: - """ + """Maps the endpoint of the unique geodesic from :math:`\mathbf{X}` to tangent vector(s) + :param grassmann_points: Point(s) on the Grassmann manifold. :param reference_point: Origin of the tangent space. :return: Point(s) on the tangent space. @@ -83,7 +85,8 @@ def log_map(grassmann_points: Union[list[Numpy2DFloatArrayOrthonormal], list[Gra @beartype def exp_map(tangent_points: list[Numpy2DFloatArray], reference_point: Union[np.ndarray, GrassmannPoint]) -> list[GrassmannPoint]: - """ + """Maps tangent vector(s) to the endpoint of a unique geodesic. + :param tangent_points: Tangent vector(s). :param reference_point: Origin of the tangent space. :return: Point(s) on the Grassmann manifold. @@ -118,7 +121,8 @@ def exp_map(tangent_points: list[Numpy2DFloatArray], def frechet_variance(grassmann_points: Union[list[Numpy2DFloatArrayOrthonormal], list[GrassmannPoint]], reference_point: Union[Numpy2DFloatArrayOrthonormal, GrassmannPoint], distance: GrassmannianDistance) -> float: - """ + """Variance measure of the Grassman distance in relation to the Karcher mean. + :param grassmann_points: Point(s) on the Grassmann manifold :param reference_point: Reference point for the Frechet variance (:math:`\mu`). Typically assigned as the Karcher mean. diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index e100cbb16..ed5b97842 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -113,6 +113,7 @@ def __init__( self.alpha: float = None """Direction cosine.""" self.failure_probability = None + """FORM probability of failure (:math:`\Phi(-\beta)`)""" self.x = None self.g0 = None self.iterations: int = None From 717328f37cae7058817e47bde085712b21066c40 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 14:47:06 -0400 Subject: [PATCH 04/22] fixing examples for updated FORM.py --- .../reliability/form/{pfn.py => local_pfn.py} | 0 .../form/plot_FORM_linear function_2d.py | 12 +++--- .../form/plot_FORM_linear_function_3d.py | 6 +-- .../form/plot_FORM_structural_reliability.py | 39 +++++++++---------- 4 files changed, 27 insertions(+), 30 deletions(-) rename docs/code/reliability/form/{pfn.py => local_pfn.py} (100%) diff --git a/docs/code/reliability/form/pfn.py b/docs/code/reliability/form/local_pfn.py similarity index 100% rename from docs/code/reliability/form/pfn.py rename to docs/code/reliability/form/local_pfn.py diff --git a/docs/code/reliability/form/plot_FORM_linear function_2d.py b/docs/code/reliability/form/plot_FORM_linear function_2d.py index 4138843ea..e9ca28ee2 100644 --- a/docs/code/reliability/form/plot_FORM_linear function_2d.py +++ b/docs/code/reliability/form/plot_FORM_linear function_2d.py @@ -17,25 +17,23 @@ #%% -import shutil - -from UQpy.run_model.RunModel import RunModel -from UQpy.run_model.model_execution.PythonModel import PythonModel from UQpy.distributions import Normal from UQpy.reliability import FORM +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel dist1 = Normal(loc=0., scale=1.) dist2 = Normal(loc=0., scale=1.) -model = PythonModel(model_script='pfn.py', model_object_name="example2") +model = PythonModel(model_script='local_pfn.py', model_object_name="example2") RunModelObject2 = RunModel(model=model) Z = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject2) Z.run() # print results -print('Design point in standard normal space: %s' % Z.DesignPoint_U) -print('Design point in original space: %s' % Z.DesignPoint_X) +print('Design point in standard normal space: %s' % Z.design_point_u) +print('Design point in original space: %s' % Z.design_point_x) print('Hasofer-Lind reliability index: %s' % Z.beta) print('FORM probability of failure: %s' % Z.failure_probability) diff --git a/docs/code/reliability/form/plot_FORM_linear_function_3d.py b/docs/code/reliability/form/plot_FORM_linear_function_3d.py index ebe3323f1..055d72728 100644 --- a/docs/code/reliability/form/plot_FORM_linear_function_3d.py +++ b/docs/code/reliability/form/plot_FORM_linear_function_3d.py @@ -34,14 +34,14 @@ dist2 = Normal(loc=5., scale=0.8) dist3 = Normal(loc=4., scale=0.4) -model = PythonModel(model_script='pfn.py', model_object_name="example3",) +model = PythonModel(model_script='local_pfn.py', model_object_name="example3",) RunModelObject3 = RunModel(model=model) Z0 = FORM(distributions=[dist1, dist2, dist3], runmodel_object=RunModelObject3) Z0.run() -print('Design point in standard normal space: %s' % Z0.DesignPoint_U) -print('Design point in original space: %s' % Z0.DesignPoint_X) +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) diff --git a/docs/code/reliability/form/plot_FORM_structural_reliability.py b/docs/code/reliability/form/plot_FORM_structural_reliability.py index 41348d524..b384f12b8 100644 --- a/docs/code/reliability/form/plot_FORM_structural_reliability.py +++ b/docs/code/reliability/form/plot_FORM_structural_reliability.py @@ -24,42 +24,41 @@ # Initially we have to import the necessary modules. #%% -import shutil -import numpy as np import matplotlib.pyplot as plt -from UQpy.run_model.RunModel import RunModel -from UQpy.run_model.model_execution.PythonModel import PythonModel +import numpy as np + from UQpy.distributions import Normal from UQpy.reliability import FORM +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel - -model = PythonModel(model_script='pfn.py', model_object_name="example1") +model = PythonModel(model_script='local_pfn.py', model_object_name="example1") RunModelObject = RunModel(model=model) dist1 = Normal(loc=200., scale=20.) -dist2 = Normal(loc=150, scale=10.) -Q = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject, tol1=1e-5, tol2=1e-5) +dist2 = Normal(loc=150., scale=10.) +Q = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject, tolerance_u=1e-5, tolerance_beta=1e-5) Q.run() - # print results -print('Design point in standard normal space: %s' % Q.DesignPoint_U) -print('Design point in original space: %s' % Q.DesignPoint_X) +print('Design point in standard normal space: %s' % Q.design_point_u) +print('Design point in original space: %s' % Q.design_point_x) print('Hasofer-Lind reliability index: %s' % Q.beta) print('FORM probability of failure: %s' % Q.failure_probability) -print(Q.dg_u_record) +print(Q.state_function_gradient_record) # Supporting function -def multivariate_gaussian(pos, mu, Sigma): +def multivariate_gaussian(pos, mu, sigma): n = mu.shape[0] - Sigma_det = np.linalg.det(Sigma) - Sigma_inv = np.linalg.inv(Sigma) - N = np.sqrt((2 * np.pi) ** n * Sigma_det) - fac = np.einsum('...k,kl,...l->...', pos - mu, Sigma_inv, pos - mu) + sigma_det = np.linalg.det(sigma) + sigma_inv = np.linalg.inv(sigma) + N = np.sqrt((2 * np.pi) ** n * sigma_det) + fac = np.einsum('...k,kl,...l->...', pos - mu, sigma_inv, pos - mu) return np.exp(-fac / 2) / N + N = 60 XX = np.linspace(150, 250, N) YX = np.linspace(120, 180, N) @@ -96,7 +95,7 @@ def multivariate_gaussian(pos, mu, Sigma): plt.rcParams.update({'font.size': 22}) plt.plot(parameters[0][0], parameters[1][0], 'r.') plt.plot([0, 200], [0, 200], 'k', linewidth=5) -plt.plot(Q.DesignPoint_X[0][0], Q.DesignPoint_X[0][1], 'bo', markersize=12) +plt.plot(Q.design_point_x[0][0], Q.design_point_x[0][1], 'bo', markersize=12) plt.contour(XX, YX, ZX, levels=20) plt.xlabel(r'$X_1$') plt.ylabel(r'$X_2$') @@ -116,9 +115,9 @@ def multivariate_gaussian(pos, mu, Sigma): plt.figure() plt.rcParams["figure.figsize"] = (12, 12) plt.rcParams.update({'font.size': 22}) -plt.plot([0, Q.DesignPoint_U[0][0]], [0, Q.DesignPoint_U[0][1]], 'b', linewidth=2) +plt.plot([0, Q.design_point_u[0][0]], [0, Q.design_point_u[0][1]], 'b', linewidth=2) plt.plot([0, -3], [5, -1], 'k', linewidth=5) -plt.plot(Q.DesignPoint_U[0][0], Q.DesignPoint_U[0][1], 'bo', markersize=12) +plt.plot(Q.design_point_u[0][0], Q.design_point_u[0][1], 'bo', markersize=12) plt.contour(XU, YU, ZU, levels=20) plt.axhline(0, color='black') plt.axvline(0, color='black') From 653558b2b09c910de158a6dfb1711669d568d74e Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 14:47:48 -0400 Subject: [PATCH 05/22] updated input and attribute names for improved FORM.py --- tests/unit_tests/reliability/test_form.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/unit_tests/reliability/test_form.py b/tests/unit_tests/reliability/test_form.py index a9bff8dd6..dadbd4576 100644 --- a/tests/unit_tests/reliability/test_form.py +++ b/tests/unit_tests/reliability/test_form.py @@ -49,7 +49,7 @@ def test_tol1_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol1=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_u=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -62,7 +62,7 @@ def test_tol2_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol2=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_beta=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -75,7 +75,7 @@ def test_tol3_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol3=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_gradient=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -88,7 +88,7 @@ def test_tol12_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol1=1.0e-3, tol2=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_u=1.0e-3, tolerance_beta=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -101,7 +101,7 @@ def test_tol13_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol1=1.0e-3, tol3=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_u=1.0e-3, tolerance_gradient=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -114,7 +114,7 @@ def test_tol23_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol3=1.0e-3, tol2=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_gradient=1.0e-3, tolerance_beta=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -127,7 +127,7 @@ def test_tol123_is_not_none(setup): dist1 = Normal(loc=200, scale=20) dist2 = Normal(loc=150, scale=10) dist = [dist1, dist2] - form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tol1=1.0e-3, tol3=1.0e-3, tol2=1.0e-3) + form_obj = FORM(distributions=dist, runmodel_object=setup, seed_u=[1, 1], tolerance_u=1.0e-3, tolerance_gradient=1.0e-3, tolerance_beta=1.0e-3) form_obj.run() for file_name in glob.glob("Model_Runs_*"): shutil.rmtree(file_name) @@ -142,13 +142,13 @@ def test_form_example(): dist1 = Normal(loc=200., scale=20.) dist2 = Normal(loc=150, scale=10.) Q = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject, - tol1=1e-5, tol2=1e-5) + tolerance_u=1e-5, tolerance_beta=1e-5) Q.run() # print results - np.allclose(Q.DesignPoint_U, np.array([-2., 1.])) - np.allclose(Q.DesignPoint_X, np.array([160., 160.])) + np.allclose(Q.design_point_u, np.array([-2., 1.])) + np.allclose(Q.design_point_x, np.array([160., 160.])) assert Q.beta[0] == 2.236067977499917 assert Q.failure_probability[0] == 0.012673659338729965 - np.allclose(Q.dg_u_record, np.array([0., 0.])) + np.allclose(Q.state_function_gradient_record, np.array([0., 0.])) From ce3d9066ebee278abbf28c4e16cf5aea3732781b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 14:49:40 -0400 Subject: [PATCH 06/22] improved documentation, cleaned up and renamed attributes, simplified variable storage after FORM.run() completes --- docs/source/reliability/form.rst | 32 +-- src/UQpy/reliability/taylor_series/FORM.py | 295 +++++++-------------- 2 files changed, 110 insertions(+), 217 deletions(-) diff --git a/docs/source/reliability/form.rst b/docs/source/reliability/form.rst index 6a4908fdd..99b3afab7 100644 --- a/docs/source/reliability/form.rst +++ b/docs/source/reliability/form.rst @@ -17,9 +17,9 @@ is the norm of the design point known as the Hasofer-Lind reliability index calc Hasofer-Lind-Rackwitz-Fiessler (HLRF) algorithm. The convergence criteria used for HLRF algorithm are: -.. math:: e1: ||\textbf{U}^{k} - \textbf{U}^{k-1}||_2 \leq 10^{-3} -.. math:: e2: ||\beta_{HL}^{k} - \beta_{HL}^{k-1}||_2 \leq 10^{-3} -.. math:: e3: ||\nabla G(\textbf{U}^{k})- \nabla G(\textbf{U}^{k-1})||_2 \leq 10^{-3} +.. math:: \text{tolerance}_U:\ ||\textbf{U}^{k} - \textbf{U}^{k-1}||_2 \leq 10^{-3} +.. math:: \text{tolerance}_\beta:\ ||\beta_{HL}^{k} - \beta_{HL}^{k-1}||_2 \leq 10^{-3} +.. math:: \text{tolerance}_{\nabla G(U)}:\ ||\nabla G(\textbf{U}^{k})- \nabla G(\textbf{U}^{k-1})||_2 \leq 10^{-3} @@ -35,29 +35,31 @@ Methods Attributes """""""""" -.. autoattribute:: UQpy.reliability.taylor_series.FORM.beta +.. autoattribute:: UQpy.reliability.taylor_series.FORM.alpha -.. autoattribute:: UQpy.reliability.taylor_series.FORM.DesignPoint_U +.. autoattribute:: UQpy.reliability.taylor_series.FORM.alpha_record -.. autoattribute:: UQpy.reliability.taylor_series.FORM.DesignPoint_X +.. autoattribute:: UQpy.reliability.taylor_series.FORM.beta -.. autoattribute:: UQpy.reliability.taylor_series.FORM.alpha +.. autoattribute:: UQpy.reliability.taylor_series.FORM.beta_record -.. autoattribute:: UQpy.reliability.taylor_series.FORM.iterations +.. autoattribute:: UQpy.reliability.taylor_series.FORM.design_point_u -.. autoattribute:: UQpy.reliability.taylor_series.FORM.u_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.design_point_x -.. autoattribute:: UQpy.reliability.taylor_series.FORM.x_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.error_record -.. autoattribute:: UQpy.reliability.taylor_series.FORM.beta_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.failure_probability -.. autoattribute:: UQpy.reliability.taylor_series.FORM.dg_u_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.iterations -.. autoattribute:: UQpy.reliability.taylor_series.FORM.alpha_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.state_function_record -.. autoattribute:: UQpy.reliability.taylor_series.FORM.g_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.state_function_gradient_record -.. autoattribute:: UQpy.reliability.taylor_series.FORM.error_record +.. autoattribute:: UQpy.reliability.taylor_series.FORM.u_record + +.. autoattribute:: UQpy.reliability.taylor_series.FORM.x_record Examples diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index ed5b97842..77f3c6f2d 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -28,46 +28,42 @@ def __init__( corr_x: Union[list, np.ndarray] = None, corr_z: Union[list, np.ndarray] = None, n_iterations: PositiveInteger = 100, - tol1: Union[float, int] = None, - tol2: Union[float, int] = None, - tol3: Union[float, int] = None, + tolerance_u: Union[float, int, None] = 1e-3, + tolerance_beta: Union[float, int, None] = 1e-3, + tolerance_gradient: Union[float, int, None] = 1e-3, ): """ - A class perform the First Order reliability Method. The :meth:`run` method of the :class:`.FORM` class can be invoked many - times and each time the results are appended to the existing ones. + A class perform the First Order reliability Method. Each time :meth:`run` is called the results are appended. This is a child class of the :class:`.TaylorSeries` class. :param distributions: Marginal probability distributions of each random variable. Must be an object of type :class:`.DistributionContinuous1D` or :class:`.JointIndependent`. :param runmodel_object: The computational model. It should be of type :class:`RunModel`. - :param seed_u: The initial starting point in the uncorrelated standard normal space **U** for the `Hasofer-Lind` - algorithm. + :param seed_x: The initial starting point in the parameter space **X** for the `Hasofer-Lind` algorithm. Either `seed_u` or `seed_x` must be provided. - Default: :code:`seed_u = (0, 0, ..., 0)` If either `seed_u` or `seed_x` is provided, then the :py:meth:`run` method will be executed automatically. Otherwise, the the :py:meth:`run` method must be executed by the user. - :param seed_x: The initial starting point in the parameter space **X** for the `Hasofer-Lind` algorithm. - Either `seed_u` or `seed_x` must be provided. + :param seed_u: The initial starting point in the uncorrelated standard normal space **U** for the `Hasofer-Lind` + algorithm. Either `seed_u` or `seed_x` must be provided. Default: :code:`seed_u = (0, 0, ..., 0)` If either `seed_u` or `seed_x` is provided, then the :py:meth:`run` method will be executed automatically. Otherwise, the the :py:meth:`run` method must be executed by the user. :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` - :param corr_z: Covariance matrix - If `corr_z` is provided, it is the correlation matrix (:math:`\mathbf{C_Z}`) of the standard normal random - vector **Z** . + :param corr_x: Covariance matrix If `corr_x` is provided, it is the correlation matrix (:math:`\mathbf{C_X}`) of the random vector **X** . :param corr_z: Covariance matrix - If `corr_x` is provided, it is the correlation matrix (:math:`\mathbf{C_X}`) of the random vector **X** . If `corr_z` is provided, it is the correlation matrix (:math:`\mathbf{C_Z}`) of the standard normal random - vector **Z** . - Default: `corr_z` is specified as the identity matrix. + vector **Z** . Default: `corr_z` is specified as the identity matrix. :param n_iterations: Maximum number of iterations for the `HLRF` algorithm. Default: :math:`100` - :param tol1: Convergence threshold for criterion `e1` of the `HLRF` algorithm. Default: :math:`1.0e-3` - :param tol2: Convergence threshold for criterion `e2` of the `HLRF` algorithm. Default: :math:`1.0e-3` - :param tol3: Convergence threshold for criterion `e3` of the `HLRF` algorithm. Default: :math:`1.0e-3` - - Any number of tolerances can be provided. Only the provided tolerances will be considered for the convergence - of the algorithm. In case none of the tolerances is provided then they are considered equal to :math:`1e-3` and - all are checked for the convergence. + :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_beta: Convergence threshold for criterion :math:`||\\beta_{HL}^{k}-\\beta_{HL}^{k-1}||_2 \leq` + :code:`tolerance_beta` of the `HLRF` algorithm. Default: :math:`1.0e-3` + :param tolerance_gradient: Convergence threshold for criterion + :math:`||\\nabla G(\mathbf{U}^{k})- \\nabla G(\mathbf{U}^{k-1})||_2 \leq` :code:`tolerance_gradient` + of the `HLRF` algorithm. Default: :math:`1.0e-3` + + By default, all three tolerances must be satisfied for convergence. + Specifying a tolerance as :code:`None` will ignore that criteria. """ if isinstance(distributions, list): self.dimension = len(distributions) @@ -92,63 +88,57 @@ def __init__( self.corr_x = corr_x self.corr_z = corr_z + self.df_step = df_step self.dist_object = distributions self.n_iterations = n_iterations self.runmodel_object = runmodel_object - self.tol1 = tol1 - self.tol2 = tol2 - self.tol3 = tol3 self.seed_u = seed_u self.seed_x = seed_x + self.tolerance_beta = tolerance_beta + self.tolerance_gradient = tolerance_gradient + self.tolerance_u = tolerance_u self.logger = logging.getLogger(__name__) # Initialize output - self.beta: float = None + self.alpha: float = np.nan + """Direction cosine.""" + self.alpha_record : list = [] + """Record of the alpha (directional cosine).""" + self.beta: list = [] """Hasofer-Lind reliability index.""" - self.DesignPoint_U: list = None + self.beta_record: list = [] + """Record of all Hasofer-Lind reliability index values.""" + self.design_point_u: list = [] """Design point in the uncorrelated standard normal space U.""" - self.DesignPoint_X: list = None + self.design_point_x: list = [] """Design point in the parameter space X.""" - self.alpha: float = None - """Direction cosine.""" - self.failure_probability = None - """FORM probability of failure (:math:`\Phi(-\beta)`)""" - self.x = None + self.error_record: list = [] + """Record of the error defined by + criteria :math:`\\text{tolerance}_U, \\text{tolerance}_\\beta, \\text{tolerance}_{\\nabla G(U)}`.""" + self.failure_probability: list = [] + """FORM probability of failure :math:`\Phi(-\\beta)`.""" self.g0 = None - self.iterations: int = None + self.iterations: list = [] """Number of model evaluations.""" - self.df_step = df_step - self.error_record: list = None - """Record of the error defined by criteria `e1, e2, e3`.""" - - self.tol1 = tol1 - self.tol2 = tol2 - self.tol3 = tol3 - - self.u_record: list = None - """Record of all iteration points in the standard normal space **U**.""" - self.x_record: list = None - """Record of all iteration points in the parameter space **X**.""" - self.g_record: list = None - """Record of the performance function.""" - self.dg_u_record: list = None - """Record of the model’s gradient in the standard normal space.""" - self.alpha_record: list = None - """Record of the alpha (directional cosine).""" - self.beta_record: list = None - """Record of all Hasofer-Lind reliability index values.""" self.jacobian_zx = None """Jacobian of the transformation from correlated standard normal space to the parameter space.""" - self.call = None + self.state_function_record: list = [] + """Record of the performance function.""" + self.state_function_gradient_record: list = [] + """Record of the model’s gradient in the standard normal space.""" + self.u_record: list = [] + """Record of all iteration points in the standard normal space **U**.""" + self.x = None + self.x_record: list = [] + """Record of all iteration points in the parameter space **X**.""" 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): + def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.ndarray] = None): """ Runs FORM. @@ -172,19 +162,14 @@ def run(self, seed_x: Union[list, np.ndarray] = None, else: raise ValueError("UQpy: Only one seed (seed_x or seed_u) must be provided") - u_record = list() - x_record = list() - g_record = list() - alpha_record = list() - error_record = list() - converged = False k = 0 beta = np.zeros(shape=(self.n_iterations + 1,)) u = np.zeros([self.n_iterations + 1, self.dimension]) + state_function_gradient_record = np.zeros([self.n_iterations + 1, self.dimension]) u[0, :] = seed - g_record.append(0.0) - dg_u_record = np.zeros([self.n_iterations + 1, self.dimension]) + + self.state_function_record.append(0.0) while not converged and k < self.n_iterations: self.logger.info("Number of iteration: %i", k) @@ -204,153 +189,59 @@ def run(self, seed_x: Union[list, np.ndarray] = None, self.jacobian_zx = self.nataf_object.jxz self.x = x - u_record.append(u) - x_record.append(x) - self.logger.info( - "Design point Y: {0}\n".format(u[k, :]) - + "Design point X: {0}\n".format(self.x) - + "Jacobian Jzx: {0}\n".format(self.jacobian_zx)) + self.u_record.append(u) + self.x_record.append(x) + self.logger.info("Design point Y: {0}\n".format(u[k, :]) + + "Design point X: {0}\n".format(self.x) + + "Jacobian Jzx: {0}\n".format(self.jacobian_zx)) # 2. evaluate Limit State Function and the gradient at point u_k and direction cosines - dg_u, qoi, _ = self._derivatives( - point_u=u[k, :], - point_x=self.x, - runmodel_object=self.runmodel_object, - nataf_object=self.nataf_object, - df_step=self.df_step, - order="first") - g_record.append(qoi) - - dg_u_record[k + 1, :] = dg_u - norm_grad = np.linalg.norm(dg_u_record[k + 1, :]) - alpha = dg_u / norm_grad - self.logger.info( - "Directional cosines (alpha): {0}\n".format(alpha) - + "Gradient (dg_y): {0}\n".format(dg_u_record[k + 1, :]) - + "norm dg_y:", - norm_grad,) - + state_function_gradient, qoi, _ = self._derivatives(point_u=u[k, :], + point_x=self.x, + runmodel_object=self.runmodel_object, + nataf_object=self.nataf_object, + df_step=self.df_step, + order="first") + self.state_function_record.append(qoi) + state_function_gradient_record[k + 1, :] = state_function_gradient + norm_of_state_function_gradient = np.linalg.norm(state_function_gradient) + alpha = state_function_gradient / norm_of_state_function_gradient + self.logger.info("Directional cosines (alpha): {0}\n".format(alpha) + + "State Function Gradient: {0}\n".format(state_function_gradient) + + "Norm of State Function Gradient: {0}\n".format(norm_of_state_function_gradient)) self.alpha = alpha.squeeze() - alpha_record.append(self.alpha) + self.alpha_record.append(self.alpha) beta[k] = -np.inner(u[k, :].T, self.alpha) - beta[k + 1] = beta[k] + qoi / norm_grad - self.logger.info( - "Beta: {0}\n".format(beta[k]) - + "Pf: {0}".format(stats.norm.cdf(-beta[k]))) + beta[k + 1] = beta[k] + qoi / norm_of_state_function_gradient + self.logger.info("Beta: {0}\n".format(beta[k]) + "Pf: {0}".format(stats.norm.cdf(-beta[k]))) u[k + 1, :] = -beta[k + 1] * self.alpha - if (self.tol1 is not None) and (self.tol2 is not None) and (self.tol3 is not None): - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error3 = np.linalg.norm(dg_u_record[k + 1, :] - dg_u_record[k, :]) - error_record.append([error1, error2, error3]) - if error1 <= self.tol1 and error2 <= self.tol2 and error3 < self.tol3: - converged = True - else: - k = k + 1 - - if (self.tol1 is None) and (self.tol2 is None) and (self.tol3 is None): - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error3 = np.linalg.norm(dg_u_record[k + 1, :] - dg_u_record[k, :]) - error_record.append([error1, error2, error3]) - if error1 <= 1e-3 or error2 <= 1e-3 or error3 < 1e-3: - converged = True - else: - k = k + 1 + error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) + error2 = np.linalg.norm(beta[k + 1] - beta[k]) + error3 = np.linalg.norm(state_function_gradient - state_function_gradient_record[k, :]) + self.error_record.append([error1, error2, error3]) - elif (self.tol1 is not None) and (self.tol2 is None) and (self.tol3 is None): - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error_record.append(error1) - if error1 <= self.tol1: - converged = True - else: - k = k + 1 - - elif ( - (self.tol1 is None) and (self.tol2 is not None) and (self.tol3 is None) - ): - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error_record.append(error2) - if error2 <= self.tol2: - converged = True - else: - k = k + 1 - - elif (self.tol1 is None) and (self.tol2 is None) and (self.tol3 is not None): - error3 = np.linalg.norm(dg_u_record[k + 1, :] - dg_u_record[k, :]) - error_record.append(error3) - if error3 < self.tol3: - converged = True - else: - k = k + 1 - - elif (self.tol1 is not None) and (self.tol2 is not None) and (self.tol3 is None): - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error_record.append([error1, error2]) - if error1 <= self.tol1 and error2 <= self.tol1: - converged = True - else: - k = k + 1 - - elif (self.tol1 is not None) and (self.tol2 is None) and (self.tol3 is not None): - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error3 = np.linalg.norm(dg_u_record[k + 1, :] - dg_u_record[k, :]) - error_record.append([error1, error3]) - if error1 <= self.tol1 and error3 < self.tol3: - converged = True - else: - k = k + 1 - - - elif (self.tol1 is None) and (self.tol2 is not None) and (self.tol3 is not None): - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error3 = np.linalg.norm(dg_u_record[k + 1, :] - dg_u_record[k, :]) - error_record.append([error2, error3]) - if error2 <= self.tol2 and error3 < self.tol3: - converged = True - else: - k = k + 1 - - self.logger.info("Error: %s", error_record[-1]) + converged_in_error1 = True if (self.tolerance_u is None) else (error1 <= self.tolerance_u) + converged_in_error2 = True if (self.tolerance_beta is None) else (error2 <= self.tolerance_beta) + converged_in_error3 = True if (self.tolerance_gradient is None) else (error3 <= self.tolerance_gradient) + converged = converged_in_error1 and converged_in_error2 and converged_in_error3 + if not converged: + k += 1 + self.logger.info("Error: %s", self.error_record[-1]) if k > self.n_iterations: self.logger.info("UQpy: Maximum number of iterations {0} was reached before convergence." .format(self.n_iterations)) - self.error_record = error_record - self.u_record = [u_record] - self.x_record = [x_record] - self.g_record = [g_record] - self.dg_u_record = [dg_u_record[:k]] - self.alpha_record = [alpha_record] else: - if self.call is None: - self.beta_record = [beta[:k]] - self.error_record = error_record - self.beta = [beta[k]] - self.DesignPoint_U = [u[k, :]] - self.DesignPoint_X = [np.squeeze(self.x)] - self.failure_probability = [stats.norm.cdf(-self.beta[-1])] - self.iterations = [k] - self.u_record = [u_record[:k]] - self.x_record = [x_record[:k]] - self.g_record = [g_record] - self.dg_u_record = [dg_u_record[:k]] - self.alpha_record = [alpha_record] - else: - self.beta_record = self.beta_record + [beta[:k]] - self.beta = self.beta + [beta[k]] - self.error_record = self.error_record + error_record - self.DesignPoint_U = self.DesignPoint_U + [u[k, :]] - self.DesignPoint_X = self.DesignPoint_X + [np.squeeze(self.x)] - self.failure_probability = self.failure_probability + [stats.norm.cdf(-beta[k])] - self.iterations = self.iterations + [k] - self.u_record = self.u_record + [u_record[:k]] - self.x_record = self.x_record + [x_record[:k]] - self.g_record = self.g_record + [g_record] - self.dg_u_record = self.dg_u_record + [dg_u_record[:k]] - self.alpha_record = self.alpha_record + [alpha_record] - self.call = True + self.design_point_u.append(u[k, :]) + self.design_point_x.append(np.squeeze(self.x)) + self.beta.append(beta[k]) + self.failure_probability.append(stats.norm.cdf(-beta[k])) + self.state_function_gradient_record.append(state_function_gradient_record[:k]) + self.iterations.append(k) + + + + From 9d2f3d0f2eadb737359572f1e00988fd0eb287db Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 16 Jun 2023 14:49:54 -0400 Subject: [PATCH 07/22] removed tabs --- docs/requirements.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index fda2af357..555c82a40 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,5 +1,5 @@ - UQpy - sphinx_autodoc_typehints - sphinx_rtd_theme - sphinx_gallery - sphinxcontrib_bibtex \ No newline at end of file +UQpy +sphinx_autodoc_typehints +sphinx_rtd_theme +sphinx_gallery +sphinxcontrib_bibtex \ No newline at end of file From c20aeba21a979211bc4d4e22a18ee834bb4c8bfd Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 21 Jun 2023 16:12:34 -0400 Subject: [PATCH 08/22] added specific versions to packages for compatability --- docs/requirements.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index 555c82a40..2346d52b8 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,5 +1,5 @@ UQpy -sphinx_autodoc_typehints -sphinx_rtd_theme -sphinx_gallery -sphinxcontrib_bibtex \ No newline at end of file +sphinx_autodoc_typehints == 1.23.0 +sphinx_rtd_theme == 1.2.0 +sphinx_gallery == 0.13.0 +sphinxcontrib_bibtex == 2.5.0 \ No newline at end of file From 1e7844b52444123e652e0675e5b8a6a1ec3ff4cd Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 21 Jun 2023 16:12:58 -0400 Subject: [PATCH 09/22] fixed minor formatting error in math statement, alphabetized class attributes --- docs/source/stochastic_process/translation.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/source/stochastic_process/translation.rst b/docs/source/stochastic_process/translation.rst index 0f18aa92f..9559b9b21 100644 --- a/docs/source/stochastic_process/translation.rst +++ b/docs/source/stochastic_process/translation.rst @@ -14,14 +14,14 @@ function of the prescribed non-Gaussian probability distribution. The nonlinear translation in the equation above has the effect of distorting the correlation of the stochastic process. That is, if the Gaussian process has correlation function :math:`\rho(\tau)` where :math:`\tau=t_2-t_1` is a time lag, -then the non-Gaussian correlation function of the process :math:`X(t)` can be compuated as: +then the non-Gaussian correlation function of the process :math:`X(t)` can be computed as: .. math:: \xi(\tau) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F^{-1}(\Phi(y_1)) F^{-1}(\Phi(y_2)) \phi(y_1, y_2; \rho(\tau)) dy_1 dy_2 where :math:`\phi(y_1, y_2; \rho(\tau))` is the bivariate normal probability density function having correlation :math:`\rho(\tau)`. This operation is known as correlation distortion and is not, in general, invertible. That is, given -the non-Gaussian correlation function :math:`\xi(\tau) ` and an arbitrarily prescribed non-Gaussian probability -distribution with cdf :math:`F(x)`, it is not always possible to identify a correponding Gaussian process having +the non-Gaussian correlation function :math:`\xi(\tau)` and an arbitrarily prescribed non-Gaussian probability +distribution with cdf :math:`F(x)`, it is not always possible to identify a corresponding Gaussian process having correlation function :math:`\rho(\tau)` that can be translated to this non-Gaussian process through the equations above :cite:`StochasticProcess11`. @@ -47,10 +47,10 @@ Methods Attributes """""""""" -.. autoattribute:: UQpy.stochastic_process.Translation.samples_gaussian .. autoattribute:: UQpy.stochastic_process.Translation.correlation_function_non_gaussian -.. autoattribute:: UQpy.stochastic_process.Translation.scaled_correlation_function_non_gaussian .. autoattribute:: UQpy.stochastic_process.Translation.power_spectrum_non_gaussian +.. autoattribute:: UQpy.stochastic_process.Translation.samples_gaussian +.. autoattribute:: UQpy.stochastic_process.Translation.scaled_correlation_function_non_gaussian Examples @@ -76,6 +76,6 @@ Methods Attributes """""""""" -.. autoattribute:: UQpy.stochastic_process.InverseTranslation.samples_gaussian +.. autoattribute:: UQpy.stochastic_process.InverseTranslation.correlation_function_gaussian .. autoattribute:: UQpy.stochastic_process.InverseTranslation.power_spectrum_gaussian -.. autoattribute:: UQpy.stochastic_process.InverseTranslation.correlation_function_gaussian \ No newline at end of file +.. autoattribute:: UQpy.stochastic_process.InverseTranslation.samples_gaussian From 3cb5fcd07a354f03a0bfa581c16c895c05ae3ffe Mon Sep 17 00:00:00 2001 From: connor-krill Date: Fri, 23 Jun 2023 10:47:55 -0400 Subject: [PATCH 10/22] Fixes issue #194 - removed duplicate import statement --- docs/code/reliability/form/local_pfn.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/code/reliability/form/local_pfn.py b/docs/code/reliability/form/local_pfn.py index db3584d02..9837b9f52 100644 --- a/docs/code/reliability/form/local_pfn.py +++ b/docs/code/reliability/form/local_pfn.py @@ -6,6 +6,7 @@ """ import numpy as np + def example1(samples=None): g = np.zeros(samples.shape[0]) for i in range(samples.shape[0]): @@ -13,10 +14,9 @@ def example1(samples=None): 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]) @@ -30,9 +30,10 @@ def example3(samples=None): for i in range(samples.shape[0]): g[i] = 6.2*samples[i, 0] - samples[i, 1]*samples[i, 2]**2 return g - + + def example4(samples=None): g = np.zeros(samples.shape[0]) for i in range(samples.shape[0]): g[i] = samples[i, 0]*samples[i, 1] - 80 - return g \ No newline at end of file + return g From 25a9f44ec2d35762e37315b53eee833b881fc01b Mon Sep 17 00:00:00 2001 From: connor-krill Date: Mon, 24 Jul 2023 17:39:51 -0400 Subject: [PATCH 11/22] fixed typo in LaTeX in docstring --- .../code/reliability/form/plot_FORM_structural_reliability.py | 4 ++-- src/UQpy/reliability/taylor_series/FORM.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/code/reliability/form/plot_FORM_structural_reliability.py b/docs/code/reliability/form/plot_FORM_structural_reliability.py index b384f12b8..a24420289 100644 --- a/docs/code/reliability/form/plot_FORM_structural_reliability.py +++ b/docs/code/reliability/form/plot_FORM_structural_reliability.py @@ -7,9 +7,9 @@ defined in a two-dimensional parameter space consisting of a resistance :math:`R` and a stress :math:`S`. The failure happens when the stress is higher than the resistance, leading to the following limit-state function: -.. math:: \textbf{X}=\{R, S\} +.. math:: \\textbf{X}=\{R, S\} -.. math:: g(\textbf{X}) = R - S +.. math:: g(\\textbf{X}) = R - S The two random variables are independent and distributed according to: diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index 77f3c6f2d..aeedfaf29 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -117,7 +117,7 @@ def __init__( """Record of the error defined by criteria :math:`\\text{tolerance}_U, \\text{tolerance}_\\beta, \\text{tolerance}_{\\nabla G(U)}`.""" self.failure_probability: list = [] - """FORM probability of failure :math:`\Phi(-\\beta)`.""" + """FORM probability of failure :math:`\\Phi(-\\beta)`.""" self.g0 = None self.iterations: list = [] """Number of model evaluations.""" From 013a6bd563b482022eb215e7a5f7f13c87423073 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Tue, 25 Jul 2023 10:16:37 +0300 Subject: [PATCH 12/22] Fixes figures in FORM structural reliability example --- .../reliability/form/plot_FORM_structural_reliability.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/docs/code/reliability/form/plot_FORM_structural_reliability.py b/docs/code/reliability/form/plot_FORM_structural_reliability.py index a24420289..00a54c824 100644 --- a/docs/code/reliability/form/plot_FORM_structural_reliability.py +++ b/docs/code/reliability/form/plot_FORM_structural_reliability.py @@ -90,7 +90,7 @@ def multivariate_gaussian(pos, mu, sigma): ZU = multivariate_gaussian(posU, mu_U, Sigma_U) # Figure 4a -plt.figure() +fig, ax = plt.subplots() plt.rcParams["figure.figsize"] = (12, 12) plt.rcParams.update({'font.size': 22}) plt.plot(parameters[0][0], parameters[1][0], 'r.') @@ -107,12 +107,12 @@ def multivariate_gaussian(pos, mu, sigma): plt.ylim([120, 200]) plt.xlim([130, 240]) plt.grid() +ax.set_aspect('equal') plt.title('Original space') -plt.axes().set_aspect('equal', 'box') plt.show() # Figure 4b -plt.figure() +fig, ax = plt.subplots() plt.rcParams["figure.figsize"] = (12, 12) plt.rcParams.update({'font.size': 22}) plt.plot([0, Q.design_point_u[0][0]], [0, Q.design_point_u[0][1]], 'b', linewidth=2) @@ -122,7 +122,6 @@ def multivariate_gaussian(pos, mu, sigma): plt.axhline(0, color='black') plt.axvline(0, color='black') plt.plot(0, 0, 'r.') - plt.xlabel(r'$U_1$') plt.ylabel(r'$U_2$') plt.text(-1.0, 1.1, '$U^\star$=({:1.2f}, {:1.2f})'.format(-2.0, 1.0), @@ -147,6 +146,6 @@ def multivariate_gaussian(pos, mu, sigma): plt.ylim([-1, 3]) plt.xlim([-3.5, 2]) plt.grid() +ax.set_aspect('equal') plt.title('Standard Normal space') -plt.axes().set_aspect('equal', 'box') plt.show() From 9319cc15acf013075f39679ac74ef89d676f6f57 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 16:57:49 +0300 Subject: [PATCH 13/22] Fixes FORM attributes inside SORM --- src/UQpy/reliability/taylor_series/SORM.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/SORM.py b/src/UQpy/reliability/taylor_series/SORM.py index dcac92934..153f285e7 100644 --- a/src/UQpy/reliability/taylor_series/SORM.py +++ b/src/UQpy/reliability/taylor_series/SORM.py @@ -104,7 +104,7 @@ def _run(self): self.dimension = self.form_object.dimension model = self.form_object.runmodel_object - dg_u_record = self.form_object.dg_u_record + dg_u_record = self.form_object.state_function_gradient_record matrix_a = np.fliplr(np.eye(self.dimension)) matrix_a[:, 0] = self.form_object.alpha @@ -126,13 +126,13 @@ def normalize(v): r1 = np.fliplr(q).T self.logger.info("UQpy: Calculating the hessian for SORM..") - hessian_g = self._derivatives(point_u=self.form_object.DesignPoint_U[-1], - point_x=self.form_object.DesignPoint_X[-1], + hessian_g = self._derivatives(point_u=self.form_object.design_point_u[-1], + point_x=self.form_object.design_point_x[-1], runmodel_object=model, nataf_object=self.form_object.nataf_object, order="second", df_step=self.df_step, - point_qoi=self.form_object.g_record[-1][-1]) + point_qoi=self.form_object.state_function_record[-1]) matrix_b = np.dot(np.dot(r1, hessian_g), r1.T) / np.linalg.norm(dg_u_record[-1]) kappa = np.linalg.eig(matrix_b[: self.dimension - 1, : self.dimension - 1]) From 79a3dd0518d6c82ddbb837b317594658c90fdf36 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 26 Jul 2023 10:10:29 -0400 Subject: [PATCH 14/22] updated local variable names to match class attribute names --- src/UQpy/reliability/taylor_series/FORM.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index aeedfaf29..a453568d5 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -103,7 +103,7 @@ def __init__( # Initialize output self.alpha: float = np.nan """Direction cosine.""" - self.alpha_record : list = [] + self.alpha_record: list = [] """Record of the alpha (directional cosine).""" self.beta: list = [] """Hasofer-Lind reliability index.""" @@ -217,15 +217,16 @@ def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.nda u[k + 1, :] = -beta[k + 1] * self.alpha - error1 = np.linalg.norm(u[k + 1, :] - u[k, :]) - error2 = np.linalg.norm(beta[k + 1] - beta[k]) - error3 = np.linalg.norm(state_function_gradient - state_function_gradient_record[k, :]) - self.error_record.append([error1, error2, error3]) + error_u = np.linalg.norm(u[k + 1, :] - u[k, :]) + error_beta = np.linalg.norm(beta[k + 1] - beta[k]) + error_gradient = np.linalg.norm(state_function_gradient - state_function_gradient_record[k, :]) + self.error_record.append([error_u, error_beta, error_gradient]) - converged_in_error1 = True if (self.tolerance_u is None) else (error1 <= self.tolerance_u) - converged_in_error2 = True if (self.tolerance_beta is None) else (error2 <= self.tolerance_beta) - converged_in_error3 = True if (self.tolerance_gradient is None) else (error3 <= self.tolerance_gradient) - converged = converged_in_error1 and converged_in_error2 and converged_in_error3 + converged_in_u = True if (self.tolerance_u is None) else (error_u <= self.tolerance_u) + converged_in_beta = True if (self.tolerance_beta is None) else (error_beta <= self.tolerance_beta) + converged_in_gradient = True if (self.tolerance_gradient is None) \ + else (error_gradient <= self.tolerance_gradient) + converged = converged_in_u and converged_in_beta and converged_in_gradient if not converged: k += 1 From c133c3d5a9dfcdf87d7ada8f860c495bdcfb8817 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 26 Jul 2023 10:10:37 -0400 Subject: [PATCH 15/22] updated local variable names to match class attribute names --- src/UQpy/reliability/taylor_series/SORM.py | 32 +++++++++++++--------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/UQpy/reliability/taylor_series/SORM.py b/src/UQpy/reliability/taylor_series/SORM.py index 153f285e7..4c01a9f14 100644 --- a/src/UQpy/reliability/taylor_series/SORM.py +++ b/src/UQpy/reliability/taylor_series/SORM.py @@ -53,9 +53,9 @@ def build_from_first_order( corr_x: Union[list, np.ndarray] = None, corr_z: Union[list, np.ndarray] = None, n_iterations: PositiveInteger = 100, - tol1: Union[float, int] = None, - tol2: Union[float, int] = None, - tol3: Union[float, int] = None, + tolerance_u: Union[float, int] = None, + tolerance_beta: Union[float, int] = None, + tolerance_gradient: Union[float, int] = None, ): """ :param distributions: Marginal probability distributions of each random variable. Must be an object of @@ -82,16 +82,20 @@ def build_from_first_order( vector **Z** . Default: `corr_z` is specified as the identity matrix. :param n_iterations: Maximum number of iterations for the `HLRF` algorithm. Default: :math:`100` - :param tol1: Convergence threshold for criterion `e1` of the `HLRF` algorithm. Default: :math:`1.0e-3` - :param tol2: Convergence threshold for criterion `e2` of the `HLRF` algorithm. Default: :math:`1.0e-3` - :param tol3: Convergence threshold for criterion `e3` of the `HLRF` algorithm. Default: :math:`1.0e-3` + :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_beta: Convergence threshold for criterion :math:`||\\beta_{HL}^{k}-\\beta_{HL}^{k-1}||_2 \leq` + :code:`tolerance_beta` of the `HLRF` algorithm. Default: :math:`1.0e-3` + :param tolerance_gradient: Convergence threshold for criterion + :math:`||\\nabla G(\mathbf{U}^{k})- \\nabla G(\mathbf{U}^{k-1})||_2 \leq` :code:`tolerance_gradient` + of the `HLRF` algorithm. Default: :math:`1.0e-3` Any number of tolerances can be provided. Only the provided tolerances will be considered for the convergence of the algorithm. In case none of the tolerances is provided then they are considered equal to :math:`1e-3` and all are checked for the convergence. """ f = FORM(distributions, runmodel_object, seed_x, seed_u, df_step, - corr_x, corr_z, n_iterations, tol1, tol2, tol3) + corr_x, corr_z, n_iterations, tolerance_u, tolerance_beta, tolerance_gradient) return cls(f, df_step) @@ -104,16 +108,17 @@ def _run(self): self.dimension = self.form_object.dimension model = self.form_object.runmodel_object - dg_u_record = self.form_object.state_function_gradient_record + state_function_gradient_record = self.form_object.state_function_gradient_record matrix_a = np.fliplr(np.eye(self.dimension)) matrix_a[:, 0] = self.form_object.alpha - def normalize(v): - return v / np.sqrt(v.dot(v)) + # def normalize(v): + # return v / np.sqrt(v.dot(v)) q = np.zeros(shape=(self.dimension, self.dimension)) - q[:, 0] = normalize(matrix_a[:, 0]) + # q[:, 0] = normalize(matrix_a[:, 0]) + q[:, 0] = matrix_a[:, 0] / np.linalg.norm() for i in range(1, self.dimension): ai = matrix_a[:, i] @@ -121,7 +126,8 @@ def normalize(v): aj = matrix_a[:, j] t = ai.dot(aj) ai = ai - t * aj - q[:, i] = normalize(ai) + # q[:, i] = normalize(ai) + q[:, i] = ai / np.linalg.norm(ai) r1 = np.fliplr(q).T self.logger.info("UQpy: Calculating the hessian for SORM..") @@ -134,7 +140,7 @@ def normalize(v): df_step=self.df_step, point_qoi=self.form_object.state_function_record[-1]) - matrix_b = np.dot(np.dot(r1, hessian_g), r1.T) / np.linalg.norm(dg_u_record[-1]) + matrix_b = np.dot(np.dot(r1, hessian_g), r1.T) / np.linalg.norm(state_function_gradient_record[-1]) kappa = np.linalg.eig(matrix_b[: self.dimension - 1, : self.dimension - 1]) if self.call is None: self.failure_probability = [stats.norm.cdf(-1 * beta) * np.prod(1 / (1 + beta * kappa[0]) ** 0.5)] From 8324ebee980e7bdf5777d63374d10d525fcaff1f Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 17:12:09 +0300 Subject: [PATCH 16/22] Fixes minor SORM error --- src/UQpy/reliability/taylor_series/SORM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/UQpy/reliability/taylor_series/SORM.py b/src/UQpy/reliability/taylor_series/SORM.py index 4c01a9f14..cc775edf3 100644 --- a/src/UQpy/reliability/taylor_series/SORM.py +++ b/src/UQpy/reliability/taylor_series/SORM.py @@ -118,7 +118,7 @@ def _run(self): q = np.zeros(shape=(self.dimension, self.dimension)) # q[:, 0] = normalize(matrix_a[:, 0]) - q[:, 0] = matrix_a[:, 0] / np.linalg.norm() + q[:, 0] = matrix_a[:, 0] / np.linalg.norm(matrix_a[:, 0]) for i in range(1, self.dimension): ai = matrix_a[:, i] From 5b618473cffd07b90dfa6749282b492c8acc54f6 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 26 Jul 2023 10:31:10 -0400 Subject: [PATCH 17/22] updated documentation to improve legibility --- docs/source/reliability/sorm.rst | 1 + src/UQpy/reliability/taylor_series/SORM.py | 35 ++++++++-------------- 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/docs/source/reliability/sorm.rst b/docs/source/reliability/sorm.rst index 8efa107f8..80f51b136 100644 --- a/docs/source/reliability/sorm.rst +++ b/docs/source/reliability/sorm.rst @@ -22,6 +22,7 @@ The :class:`.SORM` class is imported using the following command: Methods """"""" .. autoclass:: UQpy.reliability.taylor_series.SORM + :members: build_from_first_order Attributes """""""""" diff --git a/src/UQpy/reliability/taylor_series/SORM.py b/src/UQpy/reliability/taylor_series/SORM.py index 4c01a9f14..7f8ce06cc 100644 --- a/src/UQpy/reliability/taylor_series/SORM.py +++ b/src/UQpy/reliability/taylor_series/SORM.py @@ -29,12 +29,12 @@ def __init__( """ self.logger = logging.getLogger(__name__) - """Design point in the parameter space **X**.""" self.df_step = df_step - """Record of the error defined by criteria `e1, e2, e3`.""" + """Finite difference step in standard normal space.""" self.failure_probability: float = None - self.beta_second_order = None + """FORM probability of failure :math:`\\Phi(-\\beta_{HL}) \\prod_{i=1}^{n-1} (1+\\beta_{HL}\\kappa_i)^{-\\frac{1}{2}}`.""" + self.beta_second_order: list = None """Hasofer-Lind reliability index using the SORM method.""" self.call = None @@ -61,25 +61,19 @@ def build_from_first_order( :param distributions: Marginal probability distributions of each random variable. Must be an object of type :class:`.DistributionContinuous1D` or :class:`.JointIndependent`. :param runmodel_object: The computational model. It should be of type :class:`RunModel`. - :param seed_u: The initial starting point in the uncorrelated standard normal space **U** for the `Hasofer-Lind` - algorithm. + :param seed_x: The initial starting point for the `Hasofer-Lind` algorithm. Either `seed_u` or `seed_x` must be provided. - Default: :code:`seed_u = (0, 0, ..., 0)` - If either `seed_u` or `seed_x` is provided, then the :py:meth:`run` method will be executed automatically. - Otherwise, the the :py:meth:`run` method must be executed by the user. - :param seed_x: The initial starting point in the parameter space **X** for the `Hasofer-Lind` algorithm. - Either `seed_u` or `seed_x` must be provided. - If either `seed_u` or `seed_x` is provided, then the :py:meth:`run` method will be executed automatically. - Otherwise, the the :py:meth:`run` method must be executed by the user. + If `seed_u` is provided, it should be a point in the uncorrelated standard normal space of **U**. + If `seed_x` is provided, it should be a point in the parameter space of **X**. + :param seed_u: Either `seed_u` or `seed_x` must be provided. + If `seed_u` is provided, it should be a point in the uncorrelated standard normal space of **U**. + If `seed_x` is provided, it should be a point in the parameter space of **X**. :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` - :param corr_z: Covariance matrix - If `corr_z` is provided, it is the correlation matrix (:math:`\mathbf{C_Z}`) of the standard normal random - vector **Z** . + :param corr_x: Covariance matrix If `corr_x` is provided, it is the correlation matrix (:math:`\mathbf{C_X}`) of the random vector **X** . :param corr_z: Covariance matrix - If `corr_x` is provided, it is the correlation matrix (:math:`\mathbf{C_X}`) of the random vector **X** . If `corr_z` is provided, it is the correlation matrix (:math:`\mathbf{C_Z}`) of the standard normal random - vector **Z** . + vector **Z** . Default: `corr_z` is specified as the identity matrix. Default: `corr_z` is specified as the identity matrix. :param n_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` @@ -91,7 +85,7 @@ def build_from_first_order( of the `HLRF` algorithm. Default: :math:`1.0e-3` Any number of tolerances can be provided. Only the provided tolerances will be considered for the convergence - of the algorithm. In case none of the tolerances is provided then they are considered equal to :math:`1e-3` and + of the algorithm. In case none of the tolerances are provided then they are considered equal to :math:`1e-3` and all are checked for the convergence. """ f = FORM(distributions, runmodel_object, seed_x, seed_u, df_step, @@ -113,11 +107,7 @@ def _run(self): matrix_a = np.fliplr(np.eye(self.dimension)) matrix_a[:, 0] = self.form_object.alpha - # def normalize(v): - # return v / np.sqrt(v.dot(v)) - q = np.zeros(shape=(self.dimension, self.dimension)) - # q[:, 0] = normalize(matrix_a[:, 0]) q[:, 0] = matrix_a[:, 0] / np.linalg.norm() for i in range(1, self.dimension): @@ -126,7 +116,6 @@ def _run(self): aj = matrix_a[:, j] t = ai.dot(aj) ai = ai - t * aj - # q[:, i] = normalize(ai) q[:, i] = ai / np.linalg.norm(ai) r1 = np.fliplr(q).T From eb0eca1e136d0ea85928da10c316f01032113b78 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 17:42:03 +0300 Subject: [PATCH 18/22] Hides auxiliary file in SORM example --- docs/code/reliability/sorm/{pfn.py => local_pfn.py} | 0 docs/code/reliability/sorm/plot_SORM_nonlinear_function.py | 5 +---- 2 files changed, 1 insertion(+), 4 deletions(-) rename docs/code/reliability/sorm/{pfn.py => local_pfn.py} (100%) diff --git a/docs/code/reliability/sorm/pfn.py b/docs/code/reliability/sorm/local_pfn.py similarity index 100% rename from docs/code/reliability/sorm/pfn.py rename to docs/code/reliability/sorm/local_pfn.py diff --git a/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py b/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py index 6dd0bd924..255386137 100644 --- a/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py +++ b/docs/code/reliability/sorm/plot_SORM_nonlinear_function.py @@ -17,9 +17,6 @@ # Initially we have to import the necessary modules. #%% - -import shutil - import numpy as np from UQpy.run_model.RunModel import RunModel from UQpy.run_model.model_execution.PythonModel import PythonModel @@ -37,7 +34,7 @@ dist1 = Normal(loc=20., scale=2) dist2 = Lognormal(s=s, loc=0.0, scale=scale) -model = PythonModel(model_script='pfn.py', model_object_name="example4",) +model = PythonModel(model_script='local_pfn.py', model_object_name="example4",) RunModelObject4 = RunModel(model=model) form = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject4) form.run() From f15aee5dc22b2924189c4dfd9975cc57a90fbba8 Mon Sep 17 00:00:00 2001 From: connor-krill Date: Wed, 26 Jul 2023 12:05:46 -0400 Subject: [PATCH 19/22] made the FORM plots more accurate, useful, and the code more legible --- .../form/plot_FORM_structural_reliability.py | 170 ++++++++++-------- 1 file changed, 91 insertions(+), 79 deletions(-) diff --git a/docs/code/reliability/form/plot_FORM_structural_reliability.py b/docs/code/reliability/form/plot_FORM_structural_reliability.py index 00a54c824..29ff72387 100644 --- a/docs/code/reliability/form/plot_FORM_structural_reliability.py +++ b/docs/code/reliability/form/plot_FORM_structural_reliability.py @@ -19,38 +19,69 @@ .. math:: S \sim N(150, 10) """ -#%% md +# %% md # # Initially we have to import the necessary modules. -#%% +# %% -import matplotlib.pyplot as plt import numpy as np - +import matplotlib.pyplot as plt +plt.style.use('ggplot') from UQpy.distributions import Normal from UQpy.reliability import FORM 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 `local_pfn.py `_ file can be found on +# the UQpy GitHub. It contains a simple function :code:`example1` to compute the difference between the resistence and the +# stress. + +# %% + model = PythonModel(model_script='local_pfn.py', model_object_name="example1") -RunModelObject = RunModel(model=model) +runmodel_object = RunModel(model=model) + +# %% md +# +# Now we can define the resistence and stress distributions that will be passed into :code:`FORM`. +# Along with the distributions, :code:`FORM` takes in the previously defined :code:`runmodel_object` and tolerances +# for convergences. Since :code:`tolerance_gradient` is not specified in this example, it is not considered. -dist1 = Normal(loc=200., scale=20.) -dist2 = Normal(loc=150., scale=10.) -Q = FORM(distributions=[dist1, dist2], runmodel_object=RunModelObject, tolerance_u=1e-5, tolerance_beta=1e-5) -Q.run() +# %% + +distribution_resistance = Normal(loc=200., scale=20.) +distribution_stress = Normal(loc=150., scale=10.) +form = FORM(distributions=[distribution_resistance, distribution_stress], runmodel_object=runmodel_object, + tolerance_u=1e-5, tolerance_beta=1e-5) +# %% md +# +# With everything defined we are ready to run the first-order reliability method and print the results. +# The analytic solution to this problem is :math:`\bm{u}^*=(-2, 1)` with a reliability index of +# :math:`\beta_{HL}=2.2361` and a probability of failure :math:`P_{f, \text{form}} = \Phi(-\beta_{HL}) = 0.0127` -# print results -print('Design point in standard normal space: %s' % Q.design_point_u) -print('Design point in original space: %s' % Q.design_point_x) -print('Hasofer-Lind reliability index: %s' % Q.beta) -print('FORM probability of failure: %s' % Q.failure_probability) -print(Q.state_function_gradient_record) +# %% +form.run() +print('Design point in standard normal space:', form.design_point_u) +print('Design point in original space:', form.design_point_x) +print('Hasofer-Lind reliability index:', form.beta) +print('FORM probability of failure:', form.failure_probability) +print('FORM record of the function gradient:', form.state_function_gradient_record) + +# %% md +# +# This problem can be visualized in the following plots that show the FORM results in both :math:`\bm{X}` and +# :math:`\bm{U}` space. + +# %% -# Supporting function def multivariate_gaussian(pos, mu, sigma): + """Supporting function""" n = mu.shape[0] sigma_det = np.linalg.det(sigma) sigma_inv = np.linalg.inv(sigma) @@ -68,84 +99,65 @@ def multivariate_gaussian(pos, mu, sigma): YU = np.linspace(-3, 3, N) XU, YU = np.meshgrid(XU, YU) -# Mean vector and covariance matrix in the original space -parameters = [[200, 20], [150, 10]] -mu_X = np.array([parameters[0][0], parameters[1][0]]) -Sigma_X = np.array([[parameters[0][1] ** 2, 0.0], [0.0, parameters[1][1] ** 2]]) -# Mean vector and covariance matrix in the standard normal space -mu_U = np.array([0., 0.]) -Sigma_U = np.array([[1., 0.0], [0.0, 1]]) +# %% md +# +# Define the mean vector and covariance matrix in the original :math:`\textbf{X}` space and the standard normal +# :math:`\textbf{U}` space. + +# %% +mu_X = np.array([distribution_resistance.parameters['loc'], distribution_stress.parameters['loc']]) +sigma_X = np.array([[distribution_resistance.parameters['scale']**2, 0], + [0, distribution_stress.parameters['scale']**2]]) + +mu_U = np.array([0, 0]) +sigma_U = np.array([[1, 0], + [0, 1]]) # Pack X and Y into a single 3-dimensional array for the original space posX = np.empty(XX.shape + (2,)) posX[:, :, 0] = XX posX[:, :, 1] = YX -ZX = multivariate_gaussian(posX, mu_X, Sigma_X) +ZX = multivariate_gaussian(posX, mu_X, sigma_X) # Pack X and Y into a single 3-dimensional array for the standard normal space posU = np.empty(XU.shape + (2,)) posU[:, :, 0] = XU posU[:, :, 1] = YU -ZU = multivariate_gaussian(posU, mu_U, Sigma_U) +ZU = multivariate_gaussian(posU, mu_U, sigma_U) + +# %% md +# +# Plot the :code:`FORM` solution in the original :math:`\textbf{X}` space and the standard normal :math:`\text{U}` +# space. -# Figure 4a +# %% fig, ax = plt.subplots() -plt.rcParams["figure.figsize"] = (12, 12) -plt.rcParams.update({'font.size': 22}) -plt.plot(parameters[0][0], parameters[1][0], 'r.') -plt.plot([0, 200], [0, 200], 'k', linewidth=5) -plt.plot(Q.design_point_x[0][0], Q.design_point_x[0][1], 'bo', markersize=12) -plt.contour(XX, YX, ZX, levels=20) -plt.xlabel(r'$X_1$') -plt.ylabel(r'$X_2$') -plt.text(170, 182, '$X_1 - X_2=0$', - rotation=45, - horizontalalignment='center', - verticalalignment='top', - multialignment='center') -plt.ylim([120, 200]) -plt.xlim([130, 240]) -plt.grid() +ax.contour(XX, YX, ZX, + levels=20) +ax.plot([0, 200], [0, 200], + color='black', linewidth=2, label='$G(R,S)=R-S=0$', zorder=1) +ax.scatter(mu_X[0], mu_X[1], + color='black', s=64, label='Mean $(\mu_R, \mu_S)$') +ax.scatter(form.design_point_x[0][0], form.design_point_x[0][1], + color='tab:orange', marker='*', s=100, label='Design Point', zorder=2) +ax.set(xlabel='Resistence $R$', ylabel='Stress $S$', xlim=(145, 255), ylim=(115, 185)) +ax.set_title('Original $X$ Space ') ax.set_aspect('equal') -plt.title('Original space') -plt.show() +ax.legend(loc='lower right') -# Figure 4b fig, ax = plt.subplots() -plt.rcParams["figure.figsize"] = (12, 12) -plt.rcParams.update({'font.size': 22}) -plt.plot([0, Q.design_point_u[0][0]], [0, Q.design_point_u[0][1]], 'b', linewidth=2) -plt.plot([0, -3], [5, -1], 'k', linewidth=5) -plt.plot(Q.design_point_u[0][0], Q.design_point_u[0][1], 'bo', markersize=12) -plt.contour(XU, YU, ZU, levels=20) -plt.axhline(0, color='black') -plt.axvline(0, color='black') -plt.plot(0, 0, 'r.') -plt.xlabel(r'$U_1$') -plt.ylabel(r'$U_2$') -plt.text(-1.0, 1.1, '$U^\star$=({:1.2f}, {:1.2f})'.format(-2.0, 1.0), - rotation=0, - horizontalalignment='center', - verticalalignment='top', - multialignment='center') - -plt.text(-2.1, 2.05, '$20U_1 - 10U_2 + 50=0$', - rotation=63, - horizontalalignment='center', - verticalalignment='top', - multialignment='center') - -plt.text(-1.5, 0.7, r'$\overrightarrow{\beta}$', - rotation=0, - horizontalalignment='center', - verticalalignment='top', - multialignment='center') - -plt.text(0.02, -0.2, '({:1.1f}, {:1.1f})'.format(0.0, 0.0)) -plt.ylim([-1, 3]) -plt.xlim([-3.5, 2]) -plt.grid() +ax.contour(XU, YU, ZU, + levels=20, zorder=1) +ax.plot([0, -3], [5, -1], + color='black', linewidth=2, label='$G(U_1, U_2)=0$', zorder=2) +ax.arrow(0, 0, form.design_point_u[0][0], form.design_point_u[0][1], + color='tab:blue', length_includes_head=True, width=0.05, label='$\\beta=||u^*||$', zorder=2) +ax.scatter(form.design_point_u[0][0], form.design_point_u[0][1], + color='tab:orange', marker='*', s=100, label='Design Point $u^*$', zorder=2) +ax.set(xlabel='$U_1$', ylabel='$U_2$', xlim=(-3, 3), ylim=(-3, 3)) ax.set_aspect('equal') -plt.title('Standard Normal space') +ax.set_title('Standard Normal $U$ Space') +ax.legend(loc='lower right') + plt.show() From bcb171b848ccee2ec4470a7b9d71d6d421986817 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 19:40:03 +0300 Subject: [PATCH 20/22] Adds citation for FORM structural reliability example --- .../reliability/form/plot_FORM_structural_reliability.py | 8 ++++---- docs/source/bibliography.bib | 8 ++++++++ docs/source/reliability/form.rst | 6 +++--- src/UQpy/reliability/taylor_series/FORM.py | 2 +- 4 files changed, 16 insertions(+), 8 deletions(-) diff --git a/docs/code/reliability/form/plot_FORM_structural_reliability.py b/docs/code/reliability/form/plot_FORM_structural_reliability.py index 29ff72387..b52fb3da0 100644 --- a/docs/code/reliability/form/plot_FORM_structural_reliability.py +++ b/docs/code/reliability/form/plot_FORM_structural_reliability.py @@ -3,7 +3,7 @@ 3. FORM - Structural Reliability ============================================== -The benchmark problem is a simple structural reliability problem +The benchmark problem is a simple structural reliability problem (example 7.1 in :cite:`FORM_XDu`) defined in a two-dimensional parameter space consisting of a resistance :math:`R` and a stress :math:`S`. The failure happens when the stress is higher than the resistance, leading to the following limit-state function: @@ -61,7 +61,7 @@ # %% md # # With everything defined we are ready to run the first-order reliability method and print the results. -# The analytic solution to this problem is :math:`\bm{u}^*=(-2, 1)` with a reliability index of +# The analytic solution to this problem is :math:`\textbf{u}^*=(-2, 1)` with a reliability index of # :math:`\beta_{HL}=2.2361` and a probability of failure :math:`P_{f, \text{form}} = \Phi(-\beta_{HL}) = 0.0127` # %% @@ -75,8 +75,8 @@ # %% md # -# This problem can be visualized in the following plots that show the FORM results in both :math:`\bm{X}` and -# :math:`\bm{U}` space. +# This problem can be visualized in the following plots that show the FORM results in both :math:`\textbf{X}` and +# :math:`\textbf{U}` space. # %% diff --git a/docs/source/bibliography.bib b/docs/source/bibliography.bib index 39d445a81..92af2063b 100644 --- a/docs/source/bibliography.bib +++ b/docs/source/bibliography.bib @@ -879,3 +879,11 @@ @article{Kle2D author = {Zhibao Zheng and Hongzhe Dai}, keywords = {Multi-dimensional random field, Karhunen–Loève expansion, Random field simulation, Fredholm integral equation}, } + +@misc{FORM_XDu, + title = "Probabilistic Engineering Design, Chapter 7, First Order and Second Reliability Methods", + author = "Xiaoping Du", + year = 2005, + publisher = "University of Missouri – Rolla", + howpublished = "\url{https://pdesign.sitehost.iu.edu/me360/ch7.pdf}" +} diff --git a/docs/source/reliability/form.rst b/docs/source/reliability/form.rst index 99b3afab7..23bcb7c22 100644 --- a/docs/source/reliability/form.rst +++ b/docs/source/reliability/form.rst @@ -2,7 +2,7 @@ FORM ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In FORM, the performance function is linearized according to +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 @@ -17,9 +17,9 @@ is the norm of the design point known as the Hasofer-Lind reliability index calc Hasofer-Lind-Rackwitz-Fiessler (HLRF) algorithm. The convergence criteria used for HLRF algorithm are: -.. math:: \text{tolerance}_U:\ ||\textbf{U}^{k} - \textbf{U}^{k-1}||_2 \leq 10^{-3} +.. math:: \text{tolerance}_{\textbf{U}}:\ ||\textbf{U}^{k} - \textbf{U}^{k-1}||_2 \leq 10^{-3} .. math:: \text{tolerance}_\beta:\ ||\beta_{HL}^{k} - \beta_{HL}^{k-1}||_2 \leq 10^{-3} -.. math:: \text{tolerance}_{\nabla G(U)}:\ ||\nabla G(\textbf{U}^{k})- \nabla G(\textbf{U}^{k-1})||_2 \leq 10^{-3} +.. math:: \text{tolerance}_{\nabla G(\textbf{U})}:\ ||\nabla G(\textbf{U}^{k})- \nabla G(\textbf{U}^{k-1})||_2 \leq 10^{-3} diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index a453568d5..c10c895c6 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -115,7 +115,7 @@ def __init__( """Design point in the parameter space X.""" self.error_record: list = [] """Record of the error defined by - criteria :math:`\\text{tolerance}_U, \\text{tolerance}_\\beta, \\text{tolerance}_{\\nabla G(U)}`.""" + criteria :math:`\\text{tolerance}_\\textbf{U}, \\text{tolerance}_\\beta, \\text{tolerance}_{\\nabla G(\\textbf{U})}`.""" self.failure_probability: list = [] """FORM probability of failure :math:`\\Phi(-\\beta)`.""" self.g0 = None From 857502792d0ddd8785b7e78d9b649bc004263ad2 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 19:40:25 +0300 Subject: [PATCH 21/22] Corrects SORM documentation --- docs/source/reliability/sorm.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/reliability/sorm.rst b/docs/source/reliability/sorm.rst index 80f51b136..dc6b09cf9 100644 --- a/docs/source/reliability/sorm.rst +++ b/docs/source/reliability/sorm.rst @@ -4,7 +4,7 @@ SORM In SORM the performance function is approximated by a second-order Taylor series around the design point according to -.. math:: G(\textbf{U}) = G(\textbf{U}^\star) + \nabla G(\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^\intercal + \frac{1}{2}(\textbf{U}-\textbf{U}^\star)\textbf{H}(\textbf{U}-\textbf{U}^\star) +.. math:: G(\textbf{U}) = G(\textbf{U}^\star) + \nabla G(\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^\intercal + \frac{1}{2}(\textbf{U}-\textbf{U}^\star)\textbf{H}(\textbf{U}-\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^T where :math:`\textbf{H}` is the Hessian matrix of the second derivatives of :math:`G(\textbf{U})` evaluated at :math:`\textbf{U}^*`. After the design point :math:`\textbf{U}^*` is identified and the probability of failure From 83cd4b8e693bbb968a8963fc9dbcca20767ba573 Mon Sep 17 00:00:00 2001 From: Dimitris Tsapetis Date: Wed, 26 Jul 2023 19:42:59 +0300 Subject: [PATCH 22/22] Adds SORM reference --- docs/source/reliability/sorm.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/reliability/sorm.rst b/docs/source/reliability/sorm.rst index dc6b09cf9..75664e705 100644 --- a/docs/source/reliability/sorm.rst +++ b/docs/source/reliability/sorm.rst @@ -1,7 +1,7 @@ SORM ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In SORM the performance function is approximated by a second-order Taylor series around the design point according to +In SORM :cite:`FORM_XDu` the performance function is approximated by a second-order Taylor series around the design point according to .. math:: G(\textbf{U}) = G(\textbf{U}^\star) + \nabla G(\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^\intercal + \frac{1}{2}(\textbf{U}-\textbf{U}^\star)\textbf{H}(\textbf{U}-\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^T