1#! /usr/bin/env python
2
3from __future__ import print_function
4from openturns import *
5# from math import *
6
7TESTPREAMBLE()
8RandomGenerator.SetSeed(0)
9
10try:
11    #
12    # Physical model
13    #
14
15    EtatLimite = Function(
16        ['X1', 'X2'], ['G'], ["-0.5*(X1-X2)*(X1-X2) - (X1+X2)/(sqrt(2)) + 3"])
17
18    dim = EtatLimite.getInputDimension()
19    print(dim)
20
21    #
22    # probabilistic model
23    #
24
25    mean = [0.] * dim
26    sigma = [1.] * dim
27    R = IdentityMatrix(dim)
28
29    myDistribution = Normal(mean, sigma, R)
30
31    start = myDistribution.getMean()
32    Covariance = myDistribution.getCovariance()
33
34    #
35    # limit state
36    #
37
38    vect = RandomVector(myDistribution)
39
40    output = RandomVector(EtatLimite, vect)
41
42    myEvent = ThresholdEvent(output, Less(), 0.0)
43
44    #
45    # Calculation
46    #
47    #
48    # FORM/SORM Cobyla
49    myCobyla = Cobyla()
50    myCobyla.setMaximumEvaluationNumber(100 * dim)
51    myCobyla.setMaximumAbsoluteError(1.0e-10)
52    myCobyla.setMaximumRelativeError(1.0e-10)
53    myCobyla.setMaximumResidualError(1.0e-10)
54    myCobyla.setMaximumConstraintError(1.0e-10)
55
56    myAlgoC = FORM(myCobyla, myEvent, start)
57    myAlgoC2 = SORM(myCobyla, myEvent, start)
58
59    myAlgoC.run()
60    myAlgoC2.run()
61
62    resultC = FORMResult(myAlgoC.getResult())
63    resultC2 = SORMResult(myAlgoC2.getResult())
64
65    #
66    # FORM/SORM Abdo Rackwitz
67    myAbdoRackwitz = AbdoRackwitz()
68    myAbdoRackwitz.setMaximumIterationNumber(100 * dim)
69    myAbdoRackwitz.setMaximumAbsoluteError(1.0e-10)
70    myAbdoRackwitz.setMaximumRelativeError(1.0e-10)
71    myAbdoRackwitz.setMaximumResidualError(1.0e-10)
72    myAbdoRackwitz.setMaximumConstraintError(1.0e-10)
73
74    myAlgoAR = FORM(myAbdoRackwitz, myEvent, start)
75    myAlgoAR2 = SORM(myAbdoRackwitz, myEvent, start)
76
77    myAlgoAR.run()
78    # myAlgoAR2.run()
79
80    resultAR = FORMResult(myAlgoAR.getResult())
81    # resultAR2 = SORMResult(myAlgoAR2.getResult())
82
83    #
84    # Monte Carlo
85    CoV_MC = 0.5
86    myMC = MonteCarlo(myEvent)
87    myMC.setMaximumOuterSampling(100000)
88    myMC.setBlockSize(1)
89    myMC.setMaximumCoefficientOfVariation(CoV_MC)
90    myMC.run()
91
92    #
93    # LHS
94    CoV_LHS = 0.1
95    myLHS = LHS(myEvent)
96    myLHS.setMaximumOuterSampling(100000)
97    myLHS.setBlockSize(1)
98    myLHS.setMaximumCoefficientOfVariation(CoV_LHS)
99    myLHS.run()
100
101    #
102    # Outputs
103    #
104    print("")
105    print("")
106    print(
107        "************************************************************************************************")
108    print(
109        "***************************************** FORM  COBYLA *****************************************")
110    print(
111        "************************************************************************************************")
112    print("event probability = %.5e" % PfC)
113    print("generalized reliability index = %.5f" % Beta_generalizedC)
114    print(
115        "************************************************************************************************")
116    for i in range(u_starC.getDimension()):
117        print("standard space design point = %.5f" % u_starC[i])
118    print(
119        "************************************************************************************************")
120    for i in range(x_starC.getDimension()):
121        print("physical space design point = %.5f" % x_starC[i])
122    print(
123        "************************************************************************************************")
124    print("is standard point origin in failure space? ", PtC)
125    print(
126        "************************************************************************************************")
127    for i in range(gammaC.getDimension()):
128        print("importance factors = %.5f" % gammaC[i])
129    print(
130        "************************************************************************************************")
131    print("Hasofer reliability index = %.5f" % beta_hasoferC)
132    print(
133        "************************************************************************************************")
134    for i in range(SensitivityC.getSize()):
135        for j in range(SensitivityC[i].getDimension()):
136            print("Pf sensitivity = %.5f" % SensitivityC[i][j])
137    print(
138        "************************************************************************************************")
139    print("")
140    print(
141        "************************************************************************************************")
142    print(
143        "************************************** FORM ABDO RACKWITZ **************************************")
144    print(
145        "************************************************************************************************")
146    print("event probability = %.5e" % PfAR)
147    print("generalized reliability index = %.5f" % Beta_generalizedAR)
148    print(
149        "************************************************************************************************")
150    for i in range(u_starAR.getDimension()):
151        print("standard space design point = %.5f" % u_starAR[i])
152    print(
153        "************************************************************************************************")
154    for i in range(x_starAR.getDimension()):
155        print("physical space design point = %.5f" % x_starAR[i])
156    print(
157        "************************************************************************************************")
158    print("is standard point origin in failure space? ", PtAR)
159    print(
160        "************************************************************************************************")
161    for i in range(gammaAR.getDimension()):
162        print("importance factors = %.5f" % gammaAR[i])
163    print(
164        "************************************************************************************************")
165    print("Hasofer reliability index = %.5f" % beta_hasoferAR)
166    print(
167        "************************************************************************************************")
168    for i in range(SensitivityAR.getSize()):
169        for j in range(SensitivityAR[i].getDimension()):
170            print("Pf sensitivity = %.5f" % SensitivityAR[i][j])
171    print(
172        "************************************************************************************************")
173    print("")
174    print(
175        "************************************************************************************************")
176    print(
177        "***************************************** SORM  COBYLA *****************************************")
178    print(
179        "************************************************************************************************")
180    print("Breitung event probability = %.5e" % PFBreitC2)
181    print("Breitung generalized reliability index = %.5f" % BetaBreitC2)
182    print("Hohenbichler event probability = %.5e" % PFHBC2)
183    print("Hohenbichler generalized reliability index = %.5f" % BetaHBC2)
184    print("Tvedt event probability = %.5e" % PFTvedtC2)
185    print("Tvedt generalized reliability index = %.5f" % BetaTvedtC2)
186    print(
187        "************************************************************************************************")
188    for i in range(CurvC2.getDimension()):
189        print("sorted curvatures = %.5f" % cleanScalar(CurvC2[i]))
190    print(
191        "************************************************************************************************")
192    for i in range(u_starC2.getDimension()):
193        print("standard space design point = %.5f" % u_starC2[i])
194    print(
195        "************************************************************************************************")
196    for i in range(x_starC2.getDimension()):
197        print("physical space design point = %.5f" % x_starC2[i])
198    print(
199        "************************************************************************************************")
200    print(
201        "************************************************************************************************")
202    print("is standard point origin in failure space? ", PtC2)
203    print(
204        "************************************************************************************************")
205    for i in range(gammaC2.getDimension()):
206        print("importance factors = %.5f" % gammaC2[i])
207    print(
208        "************************************************************************************************")
209    print("Hasofer reliability index = %.5f" % beta_hasoferC2)
210    print(
211        "************************************************************************************************")
212    print("")
213    print(
214        "************************************************************************************************")
215    print(
216        "************************************** SORM ABDO RACKWITZ **************************************")
217    print(
218        "************************************************************************************************")
219    print("Breitung event probability = %.5e" % PFBreitAR2)
220    print("Breitung generalized reliability index = %.5f" % BetaBreitAR2)
221    print("Hohenbichler event probability = %.5e" % PFHBAR2)
222    print("Hohenbichler generalized reliability index = %.5f" % BetaHBAR2)
223    print("Tvedt event probability = %.5e" % PFTvedtAR2)
224    print("Tvedt generalized reliability index = %.5f" % BetaTvedtAR2)
225    print(
226        "************************************************************************************************")
227    for i in range(CurvAR2.getDimension()):
228        print("sorted curvatures = %.5f" % cleanScalar(CurvAR2[i]))
229    print(
230        "************************************************************************************************")
231    for i in range(u_starAR2.getDimension()):
232        print("standard space design point = %.5f" % u_starAR2[i])
233    print(
234        "************************************************************************************************")
235    for i in range(x_starAR2.getDimension()):
236        print("physical space design point = %.5f" % x_starAR2[i])
237    print(
238        "************************************************************************************************")
239    print(
240        "************************************************************************************************")
241    print("is standard point origin in failure space? ", PtAR2)
242    print(
243        "************************************************************************************************")
244    for i in range(gammaAR2.getDimension()):
245        print("importance factors = %.5f" % gammaAR2[i])
246    print(
247        "************************************************************************************************")
248    print("Hasofer reliability index = %.5f" % beta_hasoferAR2)
249    print(
250        "************************************************************************************************")
251    print("")
252    print(
253        "************************************************************************************************")
254    print(
255        "**************************************** MONTE CARLO *******************************************")
256    print(
257        "************************************************************************************************")
258    print("Pf estimation = %.5e" % PFMC)
259    print("Pf Variance estimation = %.5e" % Variance_PF_MC)
260    print("CoV = %.5f" % CVMC)
261    print("90% Confidence Interval =", "%.5e" % length90MC)
262    print("CI at 90% =[", "%.5e" % (PFMC - 0.5 * length90MC),
263          "; %.5e" % (PFMC + 0.5 * length90MC), "]")
264    print(
265        "************************************************************************************************")
266    print("")
267    print(
268        "************************************************************************************************")
269    print(
270        "******************************************* L H S **********************************************")
271    print(
272        "************************************************************************************************")
273    print("Pf estimation = %.5e" % PFLHS)
274    print("Pf Variance estimation = %.5e" % Variance_PF_LHS)
275    print("CoV = %.5f" % CVLHS)
276    print("90% Confidence Interval =", "%.5e" % length90LHS)
277    print("CI at 90% =[", "%.5e" % (PFLHS - 0.5 * length90LHS),
278          "; %.5e" % (PFLHS + 0.5 * length90LHS), "]")
279    print(
280        "************************************************************************************************")
281    print("")
282
283except:
284    import sys
285    print("t_Waarts_concave.py", sys.exc_info()[0], sys.exc_info()[1])
286