1"""
2Use the FORM - SORM algorithms
3==============================
4"""
5# %%
6# In this example we estimate a failure probability with the `FORM` algorithm on the :ref:`cantilever beam <use-case-cantilever-beam>` example. More precisely, we show how to use the associated results:
7#
8# - the design point in both physical and standard space,
9# - the probability estimation according to the FORM approximation, and the following SORM ones: Tvedt, Hohenbichler and Breitung,
10# - the Hasofer reliability index and the generalized ones evaluated from the Breitung, Tvedt and Hohenbichler approximations,
11# - the importance factors defined as the normalized director factors of the design point in the :math:`U`-space
12# - the sensitivity factors of the Hasofer reliability index and the FORM probability.
13# - the coordinates of the mean point in the standard event space.
14#
15# The coordinates of the mean point in the standard event space is:
16#
17# .. math::
18#    \frac{1}{E_1(-\beta)}\int_{\beta}^{\infty} u_1 p_1(u_1)du_1
19#
20#
21# where :math:`E_1` is the spheric univariate distribution of the standard space and :math:`\beta` is the reliability index.
22
23# %%
24# Model definition
25# ----------------
26
27# %%
28from __future__ import print_function
29from openturns.usecases import cantilever_beam as cantilever_beam
30import openturns as ot
31import openturns.viewer as viewer
32from matplotlib import pylab as plt
33ot.Log.Show(ot.Log.NONE)
34
35# %%
36# We load the model from the usecases module :
37cb = cantilever_beam.CantileverBeam()
38
39# %%
40# We use the input parameters distribution from the data class :
41distribution = cb.distribution
42distribution.setDescription(['E', 'F', 'L', 'I'])
43
44# %%
45# We define the model
46model = cb.model
47
48# %%
49# Create the event whose probability we want to estimate.
50
51# %%
52vect = ot.RandomVector(distribution)
53G = ot.CompositeRandomVector(model, vect)
54event = ot.ThresholdEvent(G, ot.Greater(), 0.3)
55event.setName("deviation")
56
57# %%
58# FORM Analysis
59# -------------
60
61# %%
62# Define a solver
63optimAlgo = ot.Cobyla()
64optimAlgo.setMaximumEvaluationNumber(1000)
65optimAlgo.setMaximumAbsoluteError(1.0e-10)
66optimAlgo.setMaximumRelativeError(1.0e-10)
67optimAlgo.setMaximumResidualError(1.0e-10)
68optimAlgo.setMaximumConstraintError(1.0e-10)
69
70# %%
71# Run FORM
72algo = ot.FORM(optimAlgo, event, distribution.getMean())
73algo.run()
74result = algo.getResult()
75
76# %%
77# Analysis of the results
78# -----------------------
79
80# %%
81# Probability
82result.getEventProbability()
83
84# %%
85# Hasofer reliability index
86result.getHasoferReliabilityIndex()
87
88# %%
89# Design point in the standard U* space.
90
91# %%
92print(result.getStandardSpaceDesignPoint())
93
94# %%
95# Design point in the physical X space.
96
97# %%
98print(result.getPhysicalSpaceDesignPoint())
99
100# %%
101# Importance factors
102graph = result.drawImportanceFactors()
103view = viewer.View(graph)
104
105# %%
106marginalSensitivity, otherSensitivity = result.drawHasoferReliabilityIndexSensitivity()
107marginalSensitivity.setLegends(["E", "F", "L", "I"])
108marginalSensitivity.setLegendPosition('bottom')
109view = viewer.View(marginalSensitivity)
110
111# %%
112marginalSensitivity, otherSensitivity = result.drawEventProbabilitySensitivity()
113marginalSensitivity.setLegends(["E", "F", "L", "I"])
114marginalSensitivity.setLegendPosition('bottom')
115view = viewer.View(marginalSensitivity)
116
117# %%
118# Error history
119optimResult = result.getOptimizationResult()
120graphErrors = optimResult.drawErrorHistory()
121graphErrors.setLegendPosition('bottom')
122graphErrors.setYMargin(0.0)
123view = viewer.View(graphErrors)
124
125# %%
126# Get additional results with SORM
127algo = ot.SORM(optimAlgo, event, distribution.getMean())
128algo.run()
129sorm_result = algo.getResult()
130
131# %%
132# Reliability index with Breitung approximation
133sorm_result.getGeneralisedReliabilityIndexBreitung()
134
135# %%
136# ... with Hohenbichler approximation
137sorm_result.getGeneralisedReliabilityIndexHohenbichler()
138
139# %%
140# .. with Tvedt approximation
141sorm_result.getGeneralisedReliabilityIndexTvedt()
142
143# %%
144# SORM probability of the event with Breitung approximation
145sorm_result.getEventProbabilityBreitung()
146
147# %%
148# ... with Hohenbichler approximation
149sorm_result.getEventProbabilityHohenbichler()
150
151# %%
152# ... with Tvedt approximation
153sorm_result.getEventProbabilityTvedt()
154
155plt.show()
156
157# %%
158# FORM analysis with finite difference gradient
159# ---------------------------------------------
160
161# %%
162# When the considered function has no analytical expression, the gradient may not be known.
163# In this case, a constant step finite difference gradient definition may be used.
164
165# %%
166
167
168def cantilever_beam_python(X):
169    E, F, L, I = X
170    return [F*L**3/(3*E*I)]
171
172
173cbPythonFunction = ot.PythonFunction(4, 1, func=cantilever_beam_python)
174epsilon = [1e-8]*4  # Here, a constant step of 1e-8 is used for every dimension
175gradStep = ot.ConstantStep(epsilon)
176cbPythonFunction.setGradient(ot.CenteredFiniteDifferenceGradient(gradStep,
177                                                                 cbPythonFunction.getEvaluation()))
178G = ot.CompositeRandomVector(cbPythonFunction, vect)
179event = ot.ThresholdEvent(G, ot.Greater(), 0.3)
180event.setName("deviation")
181
182# %%
183# However, given the different nature of the model variables, a blended (variable)
184# finite difference step may be preferable:
185# The step depends on the location in the input space
186gradStep = ot.BlendedStep(epsilon)
187cbPythonFunction.setGradient(ot.CenteredFiniteDifferenceGradient(gradStep,
188                                                                 cbPythonFunction.getEvaluation()))
189G = ot.CompositeRandomVector(cbPythonFunction, vect)
190event = ot.ThresholdEvent(G, ot.Greater(), 0.3)
191event.setName("deviation")
192
193# %%
194# We can then run the FORM analysis in the same way as before:
195algo = ot.FORM(optimAlgo, event, distribution.getMean())
196algo.run()
197result = algo.getResult()
198