1#! /usr/bin/env python
2
3from __future__ import print_function
4from openturns import *
5from math import *
6
7TESTPREAMBLE()
8RandomGenerator.SetSeed(0)
9
10
11def cleanScalar(inScalar):
12    if (fabs(inScalar) < 1.e-10):
13        inScalar = 0.0
14    return inScalar
15
16
17try:
18
19    #
20    # Physical model
21    #
22
23    dim = 26
24    inputFunction = Description(dim)
25    inputFunction[0] = "r"
26    for i in range(1, dim):
27        inputFunction[i] = "s" + str(i)
28
29    formulas = Description(1)
30    formulas[0] = "r"
31    for i in range(1, dim):
32        formulas[0] = formulas[0] + \
33            " - (" + inputFunction[i] + " ^ 2.0) / " + str(i)
34    limitState = SymbolicFunction(inputFunction, formulas)
35
36    #
37    # Probabilistic model
38    #
39
40    mean = Point(dim, 0.2)
41    mean[0] = 0.5
42
43    sigma = Point(dim, 0.1)
44
45    R = CorrelationMatrix(dim)
46    myDistribution = Normal(mean, sigma, R)
47
48    start = myDistribution.getMean()
49
50    #
51    # Limit state
52    #
53
54    vect = RandomVector(myDistribution)
55
56    output = CompositeRandomVector(limitState, vect)
57
58    myEvent = ThresholdEvent(output, Less(), 0.0)
59
60    #
61    # Calculs
62    #
63
64    #
65    # FORM/SORM Cobyla
66    myCobyla = Cobyla()
67    myCobyla.setMaximumEvaluationNumber(100 * dim)
68    myCobyla.setMaximumAbsoluteError(1.0e-10)
69    myCobyla.setMaximumRelativeError(1.0e-10)
70    myCobyla.setMaximumResidualError(1.0e-10)
71    myCobyla.setMaximumConstraintError(1.0e-10)
72
73    myAlgoC = FORM(myCobyla, myEvent, start)
74    myAlgoC2 = SORM(myCobyla, myEvent, start)
75
76    myAlgoC.run()
77    myAlgoC2.run()
78
79    resultC = FORMResult(myAlgoC.getResult())
80    resultC2 = SORMResult(myAlgoC2.getResult())
81
82    #
83    # FORM/SORM Abdo Rackwitz
84    myAbdoRackwitz = AbdoRackwitz()
85    myAbdoRackwitz.setMaximumIterationNumber(100 * dim)
86    myAbdoRackwitz.setMaximumAbsoluteError(1.0e-10)
87    myAbdoRackwitz.setMaximumRelativeError(1.0e-10)
88    myAbdoRackwitz.setMaximumResidualError(1.0e-10)
89    myAbdoRackwitz.setMaximumConstraintError(1.0e-10)
90
91    myAlgoAR = FORM(myAbdoRackwitz, myEvent, start)
92    myAlgoAR2 = SORM(myAbdoRackwitz, myEvent, start)
93
94    myAlgoAR.run()
95    myAlgoAR2.run()
96
97    resultAR = FORMResult(myAlgoAR.getResult())
98    resultAR2 = SORMResult(myAlgoAR2.getResult())
99
100    #
101    # Monte Carlo
102    CoV_MC = 0.5
103    myMC = ProbabilitySimulationAlgorithm(myEvent, MonteCarloExperiment())
104    myMC.setMaximumOuterSampling(8000)
105    myMC.setBlockSize(1)
106    myMC.setMaximumCoefficientOfVariation(CoV_MC)
107    myMC.run()
108
109    #
110    # LHS
111    CoV_LHS = 0.1
112    myLHS = LHS(myEvent)
113    myLHS.setMaximumOuterSampling(8000)
114    myLHS.setBlockSize(1)
115    myLHS.setMaximumCoefficientOfVariation(CoV_LHS)
116    myLHS.run()
117
118    #
119    #
120    # Results
121    #
122    #
123    # FORM/SORM Cobyla
124    PfC = resultC.getEventProbability()
125    Beta_generalizedC = resultC.getGeneralisedReliabilityIndex()
126    u_starC = resultC.getStandardSpaceDesignPoint()
127    x_starC = resultC.getPhysicalSpaceDesignPoint()
128    PtC = resultC.getIsStandardPointOriginInFailureSpace(
129    ) and "true" or "false"
130    gammaC = resultC.getImportanceFactors()
131    beta_hasoferC = resultC.getHasoferReliabilityIndex()
132    SensitivityC = resultC.getEventProbabilitySensitivity()
133
134    PFBreitC2 = resultC2.getEventProbabilityBreitung()
135    BetaBreitC2 = resultC2.getGeneralisedReliabilityIndexBreitung()
136    PFHBC2 = resultC2.getEventProbabilityHohenbichler()
137    BetaHBC2 = resultC2.getGeneralisedReliabilityIndexHohenbichler()
138    PFTvedtC2 = resultC2.getEventProbabilityTvedt()
139    BetaTvedtC2 = resultC2.getGeneralisedReliabilityIndexTvedt()
140    CurvC2 = resultC2.getSortedCurvatures()
141    u_starC2 = resultC2.getStandardSpaceDesignPoint()
142    x_starC2 = resultC2.getPhysicalSpaceDesignPoint()
143    PtC2 = resultC2.getIsStandardPointOriginInFailureSpace(
144    ) and "true" or "false"
145    gammaC2 = resultC2.getImportanceFactors()
146    beta_hasoferC2 = resultC2.getHasoferReliabilityIndex()
147
148    #
149    # FORM/SORM Abdo Rackwitz
150    PfAR = resultAR.getEventProbability()
151    Beta_generalizedAR = resultAR.getGeneralisedReliabilityIndex()
152    u_starAR = resultAR.getStandardSpaceDesignPoint()
153    x_starAR = resultAR.getPhysicalSpaceDesignPoint()
154    PtAR = resultAR.getIsStandardPointOriginInFailureSpace(
155    ) and "true" or "false"
156    gammaAR = resultAR.getImportanceFactors()
157    beta_hasoferAR = resultAR.getHasoferReliabilityIndex()
158    SensitivityAR = resultAR.getEventProbabilitySensitivity()
159
160    PFBreitAR2 = resultAR2.getEventProbabilityBreitung()
161    BetaBreitAR2 = resultAR2.getGeneralisedReliabilityIndexBreitung()
162    PFHBAR2 = resultAR2.getEventProbabilityHohenbichler()
163    BetaHBAR2 = resultAR2.getGeneralisedReliabilityIndexHohenbichler()
164    PFTvedtAR2 = resultAR2.getEventProbabilityTvedt()
165    BetaTvedtAR2 = resultAR2.getGeneralisedReliabilityIndexTvedt()
166    CurvAR2 = resultAR2.getSortedCurvatures()
167    u_starAR2 = resultAR2.getStandardSpaceDesignPoint()
168    x_starAR2 = resultAR2.getPhysicalSpaceDesignPoint()
169    PtAR2 = resultAR2.getIsStandardPointOriginInFailureSpace(
170    ) and "true" or "false"
171    gammaAR2 = resultAR2.getImportanceFactors()
172    beta_hasoferAR2 = resultAR2.getHasoferReliabilityIndex()
173
174    #
175    # Monte Carlo
176    ResultMC = myMC.getResult()
177    PFMC = ResultMC.getProbabilityEstimate()
178    CVMC = ResultMC.getCoefficientOfVariation()
179    Variance_PF_MC = ResultMC.getVarianceEstimate()
180    length90MC = ResultMC.getConfidenceLength(0.90)
181
182    #
183    # LHS
184    ResultLHS = myLHS.getResult()
185    PFLHS = ResultLHS.getProbabilityEstimate()
186    CVLHS = ResultLHS.getCoefficientOfVariation()
187    Variance_PF_LHS = ResultLHS.getVarianceEstimate()
188    length90LHS = ResultLHS.getConfidenceLength(0.90)
189
190    #
191
192    #
193    # Outputs
194    #
195    print("")
196    print("")
197    print(
198        "************************************************************************************************")
199    print(
200        "***************************************** FORM  COBYLA *****************************************")
201    print(
202        "************************************************************************************************")
203    print("event probability = %.5e" % PfC)
204    print("generalized reliability index = %.5f" % Beta_generalizedC)
205    print(
206        "************************************************************************************************")
207    for i in range(u_starC.getDimension()):
208        print("standard space design point = %.4f" % u_starC[i])
209    print(
210        "************************************************************************************************")
211    for i in range(x_starC.getDimension()):
212        print("physical space design point = %.5f" % x_starC[i])
213    print(
214        "************************************************************************************************")
215    print("is standard point origin in failure space? ", PtC)
216    print(
217        "************************************************************************************************")
218    for i in range(gammaC.getDimension()):
219        print("importance factors = %.5f" % gammaC[i])
220    print(
221        "************************************************************************************************")
222    print("Hasofer reliability index = %.5f" % beta_hasoferC)
223    print(
224        "************************************************************************************************")
225    for i in range(SensitivityC.getSize()):
226        for j in range(SensitivityC[i].getDimension()):
227            print("Pf sensitivity = %.5f" % SensitivityC[i][j])
228    print(
229        "************************************************************************************************")
230    print("")
231    print(
232        "************************************************************************************************")
233    print(
234        "************************************** FORM ABDO RACKWITZ **************************************")
235    print(
236        "************************************************************************************************")
237    print("event probability = %.5e" % PfAR)
238    print("generalized reliability index = %.5f" % Beta_generalizedAR)
239    print(
240        "************************************************************************************************")
241    for i in range(u_starAR.getDimension()):
242        print("standard space design point = %.4f" % u_starAR[i])
243    print(
244        "************************************************************************************************")
245    for i in range(x_starAR.getDimension()):
246        print("physical space design point = %.5f" % x_starAR[i])
247    print(
248        "************************************************************************************************")
249    print("is standard point origin in failure space? ", PtAR)
250    print(
251        "************************************************************************************************")
252    for i in range(gammaAR.getDimension()):
253        print("importance factors = %.5f" % gammaAR[i])
254    print(
255        "************************************************************************************************")
256    print("Hasofer reliability index = %.5f" % beta_hasoferAR)
257    print(
258        "************************************************************************************************")
259    for i in range(SensitivityAR.getSize()):
260        for j in range(SensitivityAR[i].getDimension()):
261            print("Pf sensitivity = %.5f" % SensitivityAR[i][j])
262    print(
263        "************************************************************************************************")
264    print("")
265    print(
266        "************************************************************************************************")
267    print(
268        "***************************************** SORM  COBYLA *****************************************")
269    print(
270        "************************************************************************************************")
271    print("Breitung event probability = %.5e" % PFBreitC2)
272    print("Breitung generalized reliability index = %.5f" % BetaBreitC2)
273    print("Hohenbichler event probability = %.5e" % PFHBC2)
274    print("Hohenbichler generalized reliability index = %.5f" % BetaHBC2)
275    print("Tvedt event probability = %.5e" % PFTvedtC2)
276    print("Tvedt generalized reliability index = %.5f" % BetaTvedtC2)
277    print(
278        "************************************************************************************************")
279    for i in range(CurvC2.getDimension()):
280        print("sorted curvatures = %.5f" % cleanScalar(CurvC2[i]))
281    print(
282        "************************************************************************************************")
283    for i in range(u_starC2.getDimension()):
284        print("standard space design point = %.4f" % u_starC2[i])
285    print(
286        "************************************************************************************************")
287    for i in range(x_starC2.getDimension()):
288        print("physical space design point = %.5f" % x_starC2[i])
289    print(
290        "************************************************************************************************")
291    print(
292        "************************************************************************************************")
293    print("is standard point origin in failure space? ", PtC2)
294    print(
295        "************************************************************************************************")
296    for i in range(gammaC2.getDimension()):
297        print("importance factors = %.5f" % gammaC2[i])
298    print(
299        "************************************************************************************************")
300    print("Hasofer reliability index = %.5f" % beta_hasoferC2)
301    print(
302        "************************************************************************************************")
303    print("")
304    print(
305        "************************************************************************************************")
306    print(
307        "************************************** SORM ABDO RACKWITZ **************************************")
308    print(
309        "************************************************************************************************")
310    print("Breitung event probability = %.5e" % PFBreitAR2)
311    print("Breitung generalized reliability index = %.5f" % BetaBreitAR2)
312    print("Hohenbichler event probability = %.5e" % PFHBAR2)
313    print("Hohenbichler generalized reliability index = %.5f" % BetaHBAR2)
314    print("Tvedt event probability = %.5e" % PFTvedtAR2)
315    print("Tvedt generalized reliability index = %.5f" % BetaTvedtAR2)
316    print(
317        "************************************************************************************************")
318    for i in range(CurvAR2.getDimension()):
319        print("sorted curvatures = %.5f" % cleanScalar(CurvAR2[i]))
320    print(
321        "************************************************************************************************")
322    for i in range(u_starAR2.getDimension()):
323        print("standard space design point = %.4f" % u_starAR2[i])
324    print(
325        "************************************************************************************************")
326    for i in range(x_starAR2.getDimension()):
327        print("physical space design point = %.5f" % x_starAR2[i])
328    print(
329        "************************************************************************************************")
330    print(
331        "************************************************************************************************")
332    print("is standard point origin in failure space? ", PtAR2)
333    print(
334        "************************************************************************************************")
335    for i in range(gammaAR2.getDimension()):
336        print("importance factors = %.5f" % gammaAR2[i])
337    print(
338        "************************************************************************************************")
339    print("Hasofer reliability index = %.5f" % beta_hasoferAR2)
340    print(
341        "************************************************************************************************")
342    print("")
343    print(
344        "************************************************************************************************")
345    print(
346        "**************************************** MONTE CARLO *******************************************")
347    print(
348        "************************************************************************************************")
349    print("Pf estimation = %.5e" % PFMC)
350    print("Pf Variance estimation = %.5e" % Variance_PF_MC)
351    print("CoV = %.5f" % CVMC)
352    print("90% Confidence Interval =", "%.5e" % length90MC)
353    print("CI at 90% =[", "%.5e" % (PFMC - 0.5 * length90MC),
354          "; %.5e" % (PFMC + 0.5 * length90MC), "]")
355    print(
356        "************************************************************************************************")
357    print("")
358    print(
359        "************************************************************************************************")
360    print(
361        "******************************************* L H S **********************************************")
362    print(
363        "************************************************************************************************")
364    print("Pf estimation = %.5e" % PFLHS)
365    print("Pf Variance estimation = %.5e" % Variance_PF_LHS)
366    print("CoV = %.5f" % CVLHS)
367    print("90% Confidence Interval =", "%.5e" % length90LHS)
368    print("CI at 90% =[", "%.5e" % (PFLHS - 0.5 * length90LHS),
369          "; %.5e" % (PFLHS + 0.5 * length90LHS), "]")
370    print(
371        "************************************************************************************************")
372    print("")
373
374except:
375    import sys
376    print("t_Waarts_25_quadratic_terms.py",
377          sys.exc_info()[0], sys.exc_info()[1])
378