Merge pull request #222 from SURGroup/Development

Development

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

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
 #

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
 #

docs/code/reliability/form/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)
 

docs/code/reliability/sorm/pfn.py → docs/code/reliability/form/local_pfn.py

@@ -6,17 +6,17 @@
 """
 import numpy as np
 
+
 def example1(samples=None):
     g = np.zeros(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])
@@ -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
+    return g

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)
 

docs/code/reliability/form/plot_FORM_structural_reliability.py

@@ -3,13 +3,13 @@
 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:
 
-.. 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:
@@ -19,47 +19,77 @@
 .. math:: S \sim N(150, 10)
 """
 
-#%% md
+# %% md
 #
 # 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
+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 <https://github.com/SURGroup/UQpy/tree/master/docs/code/reliability/sorm>`_ 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")
+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.
 
+# %%
 
-model = PythonModel(model_script='pfn.py', model_object_name="example1")
-RunModelObject = RunModel(model=model)
+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:`\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`
 
-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)
-Q.run()
+# %%
 
+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)
 
-# print results
-print('Design point in standard normal space: %s' % Q.DesignPoint_U)
-print('Design point in original space: %s' % Q.DesignPoint_X)
-print('Hasofer-Lind reliability index: %s' % Q.beta)
-print('FORM probability of failure: %s' % Q.failure_probability)
-print(Q.dg_u_record)
+# %% md
+#
+# This problem can be visualized in the following plots that show the FORM results in both :math:`\textbf{X}` and
+# :math:`\textbf{U}` space.
 
+# %%
 
-# Supporting function
-def multivariate_gaussian(pos, mu, Sigma):
+def multivariate_gaussian(pos, mu, sigma):
+    """Supporting function"""
     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)
@@ -69,85 +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)
-
-# Figure 4a
-plt.figure()
-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.DesignPoint_X[0][0], Q.DesignPoint_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()
-plt.title('Original space')
-plt.axes().set_aspect('equal', 'box')
-plt.show()
+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.
+
+# %%
+fig, ax = plt.subplots()
+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')
+ax.legend(loc='lower right')
+
+fig, ax = plt.subplots()
+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')
+ax.set_title('Standard Normal $U$ Space')
+ax.legend(loc='lower right')
 
-# Figure 4b
-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, -3], [5, -1], 'k', linewidth=5)
-plt.plot(Q.DesignPoint_U[0][0], Q.DesignPoint_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()
-plt.title('Standard Normal space')
-plt.axes().set_aspect('equal', 'box')
 plt.show()