~ubuntu-branches/ubuntu/maverick/openturns/maverick

« back to all changes in this revision

Viewing changes to doc/src/ExamplesGuide/scriptExample_beam.py

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2008-11-18 06:32:22 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20081118063222-pa0qncclrerrqkg2
Tags: 0.12.2-1
* New upstream release
* Bug fix: "New upstream release available (0.12.2)", thanks to Jerome
  Robert (Closes: #506005).
* Applied patch by J. Robert.
* debian/control: build-depends on libxml2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
 
 
2
 
 
3
from openturns import *
 
4
 
 
5
from math import *
 
6
 
 
7
from openturns_viewer import ViewImage
 
8
 
 
9
 
 
10
#######################
 
11
### Function 'deviation'
 
12
#######################
 
13
# Create here the python lines to define the implementation of the function
 
14
 
 
15
# In order to be able to use that function with the openturns library, 
 
16
# it is necessary to define a class which derives from OpenTURNSPythonFunction
 
17
 
 
18
class modelePYTHON(OpenTURNSPythonFunction) : 
 
19
  # that following method defines the input size (4) and the output size (1)
 
20
  def __init__(self) : 
 
21
     OpenTURNSPythonFunction.__init__(self,4,1)
 
22
 
 
23
  # that following method gives the implementation of modelePYTHON
 
24
  def f(self,x) : 
 
25
     E=x[0]
 
26
     F=x[1]
 
27
     L=x[2]
 
28
     I=x[3]
 
29
     return [(F*L*L*L)/(3.*E*I)]
 
30
 
 
31
# Use that function defined in the script python 
 
32
# with the openturns library
 
33
# Create a NumericalMathFunction from modelePYTHON
 
34
deviation = NumericalMathFunction(modelePYTHON())
 
35
 
 
36
 
 
37
###########################################
 
38
### Input random vector
 
39
###########################################
 
40
 
 
41
# Create the marginal distriutions of the input random vector
 
42
distributionE = Beta(0.93, 3.2, 2.8e7, 4.8e7)
 
43
distributionF = LogNormal(30000, 9000, 15000, LogNormal.MUSIGMA)
 
44
distributionL = Uniform(250, 260)
 
45
distributionI = Beta(2.5, 4.0, 3.1e2, 4.5e2)
 
46
 
 
47
# Visualize the probability density functions
 
48
 
 
49
pdfLoiE = distributionE.drawPDF()
 
50
# Change the legend
 
51
draw_E = pdfLoiE.getDrawable(0)
 
52
draw_E.setLegendName("Beta(0.93, 3.2, 2.8e7, 4.8e7)")
 
53
pdfLoiE.setDrawable(draw_E,0)
 
54
# Change the title
 
55
pdfLoiE.setTitle("PDF of E")
 
56
 
 
57
pdfLoiE.draw("distributionE_pdf", 640, 480, GraphImplementation.EPS)
 
58
Show(pdfLoiE)
 
59
 
 
60
pdfLoiF = distributionF.drawPDF()
 
61
# Change the legend
 
62
draw_F = pdfLoiF.getDrawable(0)
 
63
draw_F.setLegendName("LogNormal(30000, 9000, 15000)")
 
64
pdfLoiF.setDrawable(draw_F,0)
 
65
# Change the title
 
66
pdfLoiF.setTitle("PDF of F")
 
67
 
 
68
pdfLoiF.draw("distributionF_pdf", 640, 480, GraphImplementation.EPS)
 
69
Show(pdfLoiF)
 
70
 
 
71
pdfLoiL = distributionL.drawPDF()
 
72
# Change the legend
 
73
draw_L = pdfLoiL.getDrawable(0)
 
74
draw_L.setLegendName("Uniform(250, 260)")
 
75
pdfLoiL.setDrawable(draw_L,0)
 
76
# Change the title
 
77
pdfLoiL.setTitle("PDF of L")
 
78
 
 
79
pdfLoiL.draw("distributionL_pdf", 640, 480, GraphImplementation.EPS)
 
80
Show(pdfLoiL)
 
81
 
 
82
 
 
83
pdfLoiI = distributionI.drawPDF()
 
84
# Change the legend
 
85
draw_I = pdfLoiI.getDrawable(0)
 
86
draw_I.setLegendName("Beta(2.5, 4.0, 3.1e2, 4.5e2)")
 
87
pdfLoiI.setDrawable(draw_I,0)
 
88
# Change the title
 
89
pdfLoiI.setTitle("PDF of I")
 
90
 
 
91
pdfLoiI.draw("distributionI_pdf", 640, 480, GraphImplementation.EPS)
 
92
Show(pdfLoiI)
 
93
 
 
94
# Create the Spearman correlation matrix of the input random vector
 
95
RS = CorrelationMatrix(4)
 
96
RS[2,3] = -0.2
 
97
 
 
98
# Evaluate the correlation matrix of the Normal copula from RS
 
99
R = NormalCopula.GetNormalCorrelationFromSpearmanCorrelation(RS)
 
100
 
 
101
# Create the Normal copula parametrized by R
 
102
copuleNormal = NormalCopula(R)
 
103
 
 
104
# Create a collection of the marginal distributions
 
105
collectionMarginals = DistributionCollection(4)
 
106
collectionMarginals[0] = Distribution(distributionE)
 
107
collectionMarginals[1] = Distribution(distributionF)
 
108
collectionMarginals[2] = Distribution(distributionL)
 
109
collectionMarginals[3] = Distribution(distributionI)
 
110
 
 
111
# Create the input probability distribution of dimension 4
 
112
inputDistribution = ComposedDistribution(collectionMarginals, Copula(copuleNormal))
 
113
 
 
114
# Give a description of each component of the input distribution
 
115
inputDescription = Description(4)
 
116
inputDescription[0] = "E"
 
117
inputDescription[1] = "F"
 
118
inputDescription[2] = "L"
 
119
inputDescription[3] = "I"
 
120
inputDistribution.setDescription(inputDescription)
 
121
 
 
122
# Create the input random vector
 
123
inputRandomVector = RandomVector(inputDistribution)
 
124
 
 
125
# Create the output variable of interest
 
126
outputVariableOfInterest =  RandomVector(deviation, inputRandomVector)
 
127
 
 
128
 
 
129
##########################
 
130
### Min/Max approach Study
 
131
##########################
 
132
 
 
133
 
 
134
####################################################
 
135
# Min/Max study with deterministic experiment plane
 
136
####################################################
 
137
 
 
138
print "###################################################"
 
139
print " Min/Max study with deterministic experiment plane"
 
140
print "###################################################"
 
141
 
 
142
 
 
143
dim = deviation.getInputNumericalPointDimension()
 
144
 
 
145
# Create the structure of the experiment plane : Composite type
 
146
 
 
147
# On each direction separately, several levels are evaluated
 
148
# here,  3 levels : +/-0.5, +/-1., +/-3. from the center
 
149
levelsNumber = 3
 
150
levels = NumericalPoint(levelsNumber, 0.0)
 
151
levels[0] = 0.5
 
152
levels[1] = 1.0
 
153
levels[2] = 3.0
 
154
# Creation of the composite plane
 
155
myPlane = Composite(dim, levels)
 
156
 
 
157
# Generation of points according to the structure of the experiment plane 
 
158
# (in a reduced centered space) 
 
159
inputSample = myPlane.generate()
 
160
 
 
161
# Scaling of the structure of the experiment plane
 
162
# scaling vector for each dimension of the levels of the structure 
 
163
# to take into account the dimension of each component
 
164
# for example : the standard deviation of each component of 'inputRandomVector' 
 
165
# in case of a RandomVector
 
166
scaling = NumericalPoint(dim)
 
167
scaling[0] = sqrt(inputRandomVector.getCovariance()[0,0])
 
168
scaling[1] = sqrt(inputRandomVector.getCovariance()[1,1])
 
169
scaling[2] = sqrt(inputRandomVector.getCovariance()[2,2])
 
170
scaling[3] = sqrt(inputRandomVector.getCovariance()[3,3])
 
171
 
 
172
inputSample.scale(scaling)
 
173
 
 
174
 
 
175
# Translation of the nonReducedSample onto the center of the experiment plane
 
176
# center = mean point of the inputRandomVector distribution
 
177
center = inputRandomVector.getMean()
 
178
inputSample.translate(center)
 
179
 
 
180
# Get the number of points in the experiment plane
 
181
pointNumber = inputSample.getSize()
 
182
 
 
183
# Evaluate the ouput variable of interest on the experiment plane
 
184
outputSample = deviation(inputSample)
 
185
 
 
186
 
 
187
# Evaluate the range of the output variable of interest on that experiment plane
 
188
minValue = outputSample.getMin()
 
189
maxValue = outputSample.getMax()
 
190
 
 
191
print "From a composite  experiment plane of size = ", pointNumber
 
192
print "Levels = ", levels[0], ", ", levels[1], ", ", levels[2]
 
193
print "Min Value = ", minValue[0]
 
194
print "Max Value = ", maxValue[0]
 
195
print ""
 
196
 
 
197
###########################################################
 
198
# Min/Max study by random sampling
 
199
###########################################################
 
200
 
 
201
print "#################################"
 
202
print " Min/Max study by random sampling"
 
203
print "#################################"
 
204
 
 
205
pointNumber = 10000
 
206
print "From random sampling = ", pointNumber
 
207
outputSample2 = outputVariableOfInterest.getNumericalSample(pointNumber)
 
208
 
 
209
minValue2 = outputSample2.getMin()
 
210
maxValue2 = outputSample2.getMax()
 
211
 
 
212
print "Min Value = ", minValue2[0]
 
213
print "Max Value = ", maxValue2[0]
 
214
print ""
 
215
 
 
216
 
 
217
 
 
218
 
 
219
 
 
220
print ""
 
221
###############################################
 
222
### Random Study : central tendance of 
 
223
### the output variable of interest
 
224
###############################################
 
225
 
 
226
print "###########################################"
 
227
print "Random Study : central tendance of"
 
228
print "the output variable of interest"
 
229
print "###########################################"
 
230
print ""
 
231
 
 
232
#####################################
 
233
# Taylor variance decomposition
 
234
#####################################
 
235
 
 
236
print "##############################"
 
237
print "Taylor variance decomposition"
 
238
print "##############################"
 
239
print ""
 
240
 
 
241
# We create a quadraticCumul algorithm
 
242
myQuadraticCumul = QuadraticCumul(outputVariableOfInterest)
 
243
 
 
244
# We compute the several elements provided by the quadratic cumul algorithm
 
245
# and evaluate the number of calculus needed
 
246
nbBefr = deviation.getEvaluationCallsNumber()
 
247
 
 
248
# Mean first order
 
249
meanFirstOrder = myQuadraticCumul.getMeanFirstOrder()[0]
 
250
nbAfter1 = deviation.getEvaluationCallsNumber()
 
251
 
 
252
# Mean second order
 
253
meanSecondOrder = myQuadraticCumul.getMeanSecondOrder()[0]
 
254
nbAfter2 = deviation.getEvaluationCallsNumber()
 
255
 
 
256
# Standard deviation
 
257
stdDeviation = sqrt(myQuadraticCumul.getCovariance()[0,0]) 
 
258
nbAfter3 = deviation.getEvaluationCallsNumber()
 
259
 
 
260
print "First order mean=", myQuadraticCumul.getMeanFirstOrder()[0]
 
261
print "Evaluation calls number = ", nbAfter1 - nbBefr
 
262
print "Second order mean=", myQuadraticCumul.getMeanSecondOrder()[0]
 
263
print "Evaluation calls number = ", nbAfter2 - nbAfter1
 
264
print "Standard deviation=", sqrt(myQuadraticCumul.getCovariance()[0,0]) 
 
265
print "Evaluation calls number = ", nbAfter3 - nbAfter2
 
266
 
 
267
print  "Importance factors="
 
268
for i in range(inputRandomVector.getDimension()) :
 
269
  print inputDistribution.getDescription()[i], " = ", myQuadraticCumul.getImportanceFactors()[i]
 
270
print ""
 
271
  
 
272
 
 
273
#############################
 
274
# Random sampling
 
275
#############################
 
276
 
 
277
print "#######################"
 
278
print "Random sampling"
 
279
print "#######################"
 
280
 
 
281
size1 = 10000
 
282
output_Sample1 = outputVariableOfInterest.getNumericalSample(size1)
 
283
outputMean = output_Sample1.computeMean()
 
284
outputCovariance = output_Sample1.computeCovariance()
 
285
 
 
286
print "Sample size = ", size1
 
287
print "Mean from sample = ", outputMean[0]
 
288
print "Standard deviation from sample = ", sqrt(outputCovariance[0,0])
 
289
print ""
 
290
 
 
291
 
 
292
##########################
 
293
# Kernel Smoothing Fitting
 
294
##########################
 
295
 
 
296
 
 
297
print "##########################"
 
298
print "# Kernel Smoothing Fitting"
 
299
print "##########################"
 
300
 
 
301
# We generate a sample of the output variable
 
302
size = 10000
 
303
output_sample = outputVariableOfInterest.getNumericalSample(size)
 
304
 
 
305
# We build the kernel smoothing distribution
 
306
kernel = KernelSmoothing()
 
307
smoothed = kernel.buildImplementation(output_sample)
 
308
print "Sample size = ", size
 
309
print  "Kernel bandwidth=" , kernel.getBandwidth()[0]
 
310
 
 
311
# We draw the pdf and cdf from kernel smoothing
 
312
# Evaluate at best the range of the graph
 
313
mean_sample = output_sample.computeMean()[0]
 
314
standardDeviation_sample = sqrt(output_sample.computeCovariance()[0,0])
 
315
xmin = mean_sample - 4*standardDeviation_sample
 
316
xmax = mean_sample + 4*standardDeviation_sample
 
317
 
 
318
# Draw the PDF
 
319
smoothedPDF = smoothed.drawPDF(xmin, xmax, 251)
 
320
# Change the title
 
321
smoothedPDF.setTitle("Kernel smoothing of the deviation - PDF")
 
322
# Change the legend
 
323
smoothedPDF_draw = smoothedPDF.getDrawable(0)
 
324
title = "PDF from Normal kernel (" + str(size) + " data)"
 
325
smoothedPDF_draw.setLegendName(title)
 
326
smoothedPDF.setDrawable(smoothedPDF_draw,0)
 
327
smoothedPDF.draw("smoothedPDF", 640, 480, GraphImplementation.EPS)
 
328
 
 
329
# Draw the CDF
 
330
smoothedCDF = smoothed.drawCDF(xmin, xmax, 251)
 
331
# Change the title
 
332
smoothedCDF.setTitle("Kernel smoothing of the deviation - CDF")
 
333
# Change the legend
 
334
smoothedCDF_draw = smoothedCDF.getDrawable(0)
 
335
title = "CDF from Normal kernel (" + str(size) + " data)"
 
336
smoothedCDF_draw.setLegendName(title)
 
337
smoothedCDF.setDrawable(smoothedCDF_draw,0)
 
338
# Change the legend position
 
339
smoothedCDF.setLegendPosition("bottomright")
 
340
smoothedCDF.draw("smoothedCDF", 640, 480, GraphImplementation.EPS)
 
341
 
 
342
# In order to see the graph whithout creating the associated files
 
343
Show(smoothedCDF) 
 
344
Show(smoothedPDF) 
 
345
 
 
346
# Mean of the output variable of interest
 
347
print "Mean from kernel smoothing = ", smoothed.getMean()[0]
 
348
print ""
 
349
 
 
350
# Superposition of the kernel smoothing pdf and the gaussian one 
 
351
# which mean and standard deviation are those of the output_sample
 
352
normalDist = NormalFactory().buildImplementation(output_sample)
 
353
normalDistPDF = normalDist.drawPDF(xmin, xmax, 251)
 
354
normalDistPDFDrawable = normalDistPDF.getDrawable(0)
 
355
normalDistPDFDrawable.setColor('blue')
 
356
smoothedPDF.addDrawable(normalDistPDFDrawable)
 
357
smoothedPDF.draw("smoothedPDF_and_NormalPDF", 640, 480, GraphImplementation.EPS)
 
358
 
 
359
# In order to see the graph whithout creating the associated files
 
360
Show(smoothedPDF) 
 
361
 
 
362
#################################################################
 
363
### Probabilistic Study : threshold exceedance: deviation > 30cm
 
364
#################################################################
 
365
 
 
366
print ""
 
367
print "############################################################"
 
368
print "Probabilistic Study : threshold exceedance: deviation <-1cm"
 
369
print "############################################################"
 
370
print ""
 
371
 
 
372
######
 
373
# FORM
 
374
######
 
375
 
 
376
print "#####"
 
377
print "FORM"
 
378
print "#####"
 
379
 
 
380
# We create an Event from this RandomVector
 
381
# threshold has been defined in the kernel smoothing section
 
382
threshold = 30
 
383
myEvent = Event(outputVariableOfInterest, ComparisonOperator(Greater()), threshold)
 
384
myEvent.setName("Deviation > 30 cm")
 
385
 
 
386
# We create a NearestPoint algorithm 
 
387
myCobyla = Cobyla()
 
388
myCobyla.setMaximumIterationsNumber(1000)
 
389
myCobyla.setMaximumAbsoluteError(1.0e-10)
 
390
myCobyla.setMaximumRelativeError(1.0e-10)
 
391
myCobyla.setMaximumResidualError(1.0e-10)
 
392
myCobyla.setMaximumConstraintError(1.0e-10)
 
393
 
 
394
# We create a FORM algorithm 
 
395
# The first parameter is a NearestPointAlgorithm 
 
396
# The second parameter is an event 
 
397
# The third parameter is a starting point for the design point research
 
398
meanVector = inputRandomVector.getMean()
 
399
myAlgoFORM = FORM(NearestPointAlgorithm(myCobyla), myEvent, meanVector)
 
400
 
 
401
# Get the number of times the limit state function has been evaluated so far
 
402
deviationCallNumberBeforeFORM = deviation.getEvaluationCallsNumber()
 
403
 
 
404
# Perform the simulation 
 
405
myAlgoFORM.run()
 
406
 
 
407
# Get the number of times the limit state function has been evaluated so far
 
408
deviationCallNumberAfterFORM = deviation.getEvaluationCallsNumber()
 
409
 
 
410
# Stream out the result 
 
411
resultFORM = myAlgoFORM.getResult()
 
412
print  "FORM event probability=" , resultFORM.getEventProbability()
 
413
print "Number of evaluations of the limit state function = ", deviationCallNumberAfterFORM - deviationCallNumberBeforeFORM
 
414
print  "Generalized reliability index=" , resultFORM.getGeneralisedReliabilityIndex() 
 
415
print  "Standard space design point="
 
416
for i in range(inputRandomVector.getDimension()) :
 
417
  print inputDistribution.getDescription()[i], " = ", resultFORM.getStandardSpaceDesignPoint()[i]
 
418
print  "Physical space design point="
 
419
for i in range(inputRandomVector.getDimension()) :
 
420
  print inputDistribution.getDescription()[i], " = ", resultFORM.getPhysicalSpaceDesignPoint()[i]
 
421
 
 
422
print  "Importance factors="
 
423
for i in range(inputRandomVector.getDimension()) :
 
424
  print inputDistribution.getDescription()[i], " = ", resultFORM.getImportanceFactors()[i]
 
425
  
 
426
print  "Hasofer reliability index=" , resultFORM.getHasoferReliabilityIndex() 
 
427
print ""
 
428
 
 
429
# Graph 1 : Importance Factors graph */
 
430
importanceFactorsGraph = resultFORM.drawImportanceFactors()
 
431
title = "FORM Importance factors - "+ myEvent.getName()
 
432
importanceFactorsGraph.setTitle( title)
 
433
importanceFactorsGraph.draw("ImportanceFactorsDrawingFORM", 640, 480, GraphImplementation.EPS)
 
434
    
 
435
# In order to see the graph whithout creating the associated files
 
436
Show(importanceFactorsGraph) 
 
437
 
 
438
 
 
439
######
 
440
# MC
 
441
######
 
442
 
 
443
print "############"
 
444
print "Monte Carlo"
 
445
print "############"
 
446
print ""
 
447
 
 
448
 
 
449
maximumOuterSampling = 40000
 
450
blockSize = 100
 
451
coefficientOfVariation = 0.10
 
452
 
 
453
# We create a Monte Carlo algorithm 
 
454
myAlgoMonteCarlo = MonteCarlo(myEvent)
 
455
myAlgoMonteCarlo.setMaximumOuterSampling(maximumOuterSampling)
 
456
myAlgoMonteCarlo.setBlockSize(blockSize)
 
457
myAlgoMonteCarlo.setMaximumCoefficientOfVariation(coefficientOfVariation)
 
458
 
 
459
# Define the HistoryStrategy to store the numerical samples generated
 
460
# both for the input random vector and the output random vector
 
461
# Full strategy 
 
462
myAlgoMonteCarlo.setInputOutputStrategy(HistoryStrategy(Full()))
 
463
 
 
464
# Define the HistoryStrategy to store the values of the probability estimator
 
465
# and the variance estimator
 
466
# used ot draw the convergence graph
 
467
# Full strategy 
 
468
myAlgoMonteCarlo.setConvergenceStrategy(HistoryStrategy(Full()))
 
469
                                  
 
470
# Perform the simulation 
 
471
myAlgoMonteCarlo.run()
 
472
 
 
473
# Display number of iterations and number of evaluations 
 
474
# of the limit state function
 
475
print "Number of evaluations of the limit state function = ", myAlgoMonteCarlo.getResult().getOuterSampling()* myAlgoMonteCarlo.getResult().getBlockSize()
 
476
 
 
477
# Display the Monte Carlo probability of 'myEvent'
 
478
print "Monte Carlo probability estimation = ", myAlgoMonteCarlo.getResult().getProbabilityEstimate()
 
479
 
 
480
# Display the variance of the Monte Carlo probability estimator
 
481
print "Variance of the Monte Carlo probability estimator = ", myAlgoMonteCarlo.getResult().getVarianceEstimate()
 
482
 
 
483
# Display the confidence interval length centered around
 
484
# the MonteCarlo probability MCProb
 
485
# IC = [MCProb - 0.5*length, MCProb + 0.5*length]
 
486
# level 0.95
 
487
 
 
488
print "0.95 Confidence Interval = [", myAlgoMonteCarlo.getResult().getProbabilityEstimate() - 0.5*myAlgoMonteCarlo.getResult().getConfidenceLength(0.95), ", ", myAlgoMonteCarlo.getResult().getProbabilityEstimate() + 0.5*myAlgoMonteCarlo.getResult().getConfidenceLength(0.95), "]"
 
489
print ""
 
490
 
 
491
# Draw the convergence graph and the confidence intervalle of level alpha
 
492
alpha = 0.90
 
493
convergenceGraphMonteCarlo = myAlgoMonteCarlo.drawProbabilityConvergence(alpha)
 
494
# In order to see the graph whithout creating the associated files
 
495
Show(convergenceGraphMonteCarlo) 
 
496
 
 
497
# Create the file .EPS
 
498
convergenceGraphMonteCarlo.draw("convergenceGrapheMonteCarlo", 640, 480, GraphImplementation.EPS)
 
499
Show(convergenceGraphMonteCarlo)
 
500
 
 
501
 
 
502
########################
 
503
# Directional Sampling
 
504
########################
 
505
 
 
506
print "#######################"
 
507
print "Directional Sampling"
 
508
print "#######################"
 
509
 
 
510
# Directional sampling from an event (slow and safe strategy by default)
 
511
 
 
512
# We create a Directional Sampling algorithm */
 
513
myAlgoDirectionalSim = DirectionalSampling(myEvent)
 
514
myAlgoDirectionalSim.setMaximumOuterSampling(maximumOuterSampling * blockSize)
 
515
myAlgoDirectionalSim.setBlockSize(1)
 
516
myAlgoDirectionalSim.setMaximumCoefficientOfVariation(coefficientOfVariation)
 
517
 
 
518
# Define the HistoryStrategy to store the numerical samples generated
 
519
# both for the input random vector and the output random vector
 
520
# Full strategy 
 
521
myAlgoDirectionalSim.setInputOutputStrategy(HistoryStrategy(Full()))
 
522
 
 
523
# Define the HistoryStrategy to store the values of the probability estimator
 
524
# and the variance estimator
 
525
# used ot draw the convergence graph
 
526
# Full strategy 
 
527
myAlgoDirectionalSim.setConvergenceStrategy(HistoryStrategy(Full()))
 
528
 
 
529
# Save the number of calls to the limit state fucntion, its gradient and hessian already done
 
530
deviationCallNumberBefore = deviation.getEvaluationCallsNumber()
 
531
deviationGradientCallNumberBefore = deviation.getGradientCallsNumber()
 
532
deviationHessianCallNumberBefore = deviation.getHessianCallsNumber()
 
533
 
 
534
# Perform the simulation */
 
535
myAlgoDirectionalSim.run()
 
536
    
 
537
# Save the number of calls to the limit state fucntion, its gradient and hessian already done
 
538
deviationCallNumberAfter = deviation.getEvaluationCallsNumber()
 
539
deviationGradientCallNumberAfter = deviation.getGradientCallsNumber()
 
540
deviationHessianCallNumberAfter = deviation.getHessianCallsNumber()
 
541
 
 
542
# Display number of iterations and number of evaluations 
 
543
# of the limit state function
 
544
print "Number of evaluations of the limit state function = ", deviationCallNumberAfter - deviationCallNumberBefore
 
545
 
 
546
# Display the Directional Simumation probability of 'myEvent'
 
547
print "Directional Sampling probability estimation = ", myAlgoDirectionalSim.getResult().getProbabilityEstimate()
 
548
 
 
549
# Display the variance of the Directional Simumation probability estimator
 
550
print "Variance of the Directional Sampling probability estimator = ", myAlgoDirectionalSim.getResult().getVarianceEstimate()
 
551
 
 
552
# Display the confidence interval length centered around 
 
553
# the Directional Simumation probability DSProb
 
554
# IC = [DSProb - 0.5*length, DSProb + 0.5*length]
 
555
# level 0.95
 
556
print "0.95 Confidence Interval = [", myAlgoDirectionalSim.getResult().getProbabilityEstimate() - 0.5*myAlgoDirectionalSim.getResult().getConfidenceLength(0.95), ", ", myAlgoDirectionalSim.getResult().getProbabilityEstimate() + 0.5*myAlgoDirectionalSim.getResult().getConfidenceLength(0.95), "]"
 
557
print ""
 
558
 
 
559
 
 
560
# Draw the convergence graph and the confidence intervalle of level alpha
 
561
alpha = 0.90
 
562
convergenceGraphDS = myAlgoDirectionalSim.drawProbabilityConvergence(alpha)
 
563
# In order to see the graph whithout creating the associated files
 
564
Show(convergenceGraphDS) 
 
565
 
 
566
# Create the file .EPS
 
567
convergenceGraphDS.draw("convergenceGrapheDS", 640, 480, GraphImplementation.EPS)
 
568
Show(convergenceGraphDS)
 
569
 
 
570
 
 
571
##########################
 
572
# Latin HyperCube Sampling
 
573
###########################
 
574
 
 
575
print "###########################"
 
576
print "Latin HyperCube Sampling"
 
577
print "###########################"
 
578
 
 
579
# We create a LHS algorithm 
 
580
myAlgoLHS = LHS(myEvent)
 
581
myAlgoLHS.setMaximumOuterSampling(maximumOuterSampling)
 
582
myAlgoLHS.setBlockSize(blockSize)
 
583
myAlgoLHS.setMaximumCoefficientOfVariation(coefficientOfVariation)
 
584
 
 
585
# Define the HistoryStrategy to store the numerical samples generated
 
586
# both for the input random vector and the output random vector
 
587
# Full strategy 
 
588
myAlgoLHS.setInputOutputStrategy(HistoryStrategy(Full()))
 
589
 
 
590
# Define the HistoryStrategy to store the values of the probability estimator
 
591
# and the variance estimator
 
592
# used ot draw the convergence graph
 
593
# Full strategy 
 
594
myAlgoLHS.setConvergenceStrategy(HistoryStrategy(Full()))
 
595
 
 
596
# Perform the simulation 
 
597
myAlgoLHS.run()
 
598
 
 
599
# Display number of iterations and number of evaluations 
 
600
# of the limit state function
 
601
print "Number of evaluations of the limit state function = ", myAlgoLHS.getResult().getOuterSampling()*myAlgoLHS.getResult().getBlockSize()
 
602
 
 
603
# Display the LHS probability of {\itshape myEvent}
 
604
print "LHS probability estimation = ", myAlgoLHS.getResult().getProbabilityEstimate()
 
605
 
 
606
# Display the variance of the LHS probability estimator
 
607
print "Variance of the LHS probability estimator = ", myAlgoLHS.getResult().getVarianceEstimate()
 
608
 
 
609
# Display the confidence interval length centered aroung the LHS probability LHSProb
 
610
# IC = [LHSProb - 0.5*length, LHSProb + 0.5*length]
 
611
# level 0.95
 
612
print "0.95 Confidence Interval = [", myAlgoLHS.getResult().getProbabilityEstimate() - 0.5*myAlgoLHS.getResult().getConfidenceLength(0.95), ", ", myAlgoLHS.getResult().getProbabilityEstimate() + 0.5*myAlgoLHS.getResult().getConfidenceLength(0.95), "]"
 
613
print ""
 
614
 
 
615
# Draw the convergence graph and the confidence intervalle of level alpha
 
616
alpha = 0.90
 
617
convergenceGraphLHS = myAlgoLHS.drawProbabilityConvergence(alpha)
 
618
# In order to see the graph whithout creating the associated files
 
619
Show(convergenceGraphLHS) 
 
620
 
 
621
# Create the file .EPS
 
622
convergenceGraphLHS.draw("convergenceGrapheLHS", 640, 480, GraphImplementation.EPS)
 
623
Show(convergenceGraphLHS)
 
624
 
 
625
#####################
 
626
# Importance Sampling
 
627
#####################
 
628
 
 
629
 
 
630
print "####################"
 
631
print "Importance Sampling"
 
632
print "####################"
 
633
 
 
634
maximumOuterSampling = 40000
 
635
blockSize = 1
 
636
standardSpaceDesignPoint = resultFORM.getStandardSpaceDesignPoint()
 
637
mean = standardSpaceDesignPoint
 
638
sigma = NumericalPoint(4, 1.0)
 
639
importanceDistribution = Normal(mean, sigma, CorrelationMatrix(4))
 
640
 
 
641
myStandardEvent = StandardEvent(myEvent)
 
642
 
 
643
myAlgoImportanceSampling = ImportanceSampling(myStandardEvent, Distribution(importanceDistribution))
 
644
myAlgoImportanceSampling.setMaximumOuterSampling(maximumOuterSampling)
 
645
myAlgoImportanceSampling.setBlockSize(blockSize)
 
646
myAlgoImportanceSampling.setMaximumCoefficientOfVariation(coefficientOfVariation)
 
647
 
 
648
# Define the HistoryStrategy to store the numerical samples generated
 
649
# both for the input random vector and the output random vector
 
650
# Full strategy 
 
651
myAlgoImportanceSampling.setInputOutputStrategy(HistoryStrategy(Full()))
 
652
 
 
653
# Define the HistoryStrategy to store the values of the probability estimator
 
654
# and the variance estimator
 
655
# used ot draw the convergence graph
 
656
# Full strategy 
 
657
myAlgoImportanceSampling.setConvergenceStrategy(HistoryStrategy(Full()))
 
658
 
 
659
# Perform the simulation 
 
660
myAlgoImportanceSampling.run()
 
661
 
 
662
# Display number of iterations and number of evaluations 
 
663
# of the limit state function
 
664
print "Number of evaluations of the limit state function = ", myAlgoImportanceSampling.getResult().getOuterSampling()* myAlgoImportanceSampling.getResult().getBlockSize()
 
665
 
 
666
# Display the Importance Sampling probability of 'myEvent'
 
667
print "Importance Sampling probability estimation = ", myAlgoImportanceSampling.getResult().getProbabilityEstimate()
 
668
 
 
669
# Display the variance of the Importance Sampling probability estimator
 
670
print "Variance of the Importance Sampling probability estimator = ", myAlgoImportanceSampling.getResult().getVarianceEstimate()
 
671
 
 
672
# Display the confidence interval length centered around
 
673
# the ImportanceSampling probability ISProb
 
674
# IC = [ISProb - 0.5*length, ISProb + 0.5*length]
 
675
# level 0.95
 
676
print "0.95 Confidence Interval = [", myAlgoImportanceSampling.getResult().getProbabilityEstimate() - 0.5*myAlgoImportanceSampling.getResult().getConfidenceLength(0.95), ", ", myAlgoImportanceSampling.getResult().getProbabilityEstimate() + 0.5*myAlgoImportanceSampling.getResult().getConfidenceLength(0.95), "]"
 
677
 
 
678
# Draw the convergence graph and the confidence intervalle of level alpha
 
679
alpha = 0.90
 
680
convergenceGraphIS = myAlgoImportanceSampling.drawProbabilityConvergence(alpha)
 
681
# In order to see the graph whithout creating the associated files
 
682
Show(convergenceGraphIS) 
 
683
 
 
684
# Create the file .EPS
 
685
convergenceGraphIS.draw("convergenceGrapheIS", 640, 480, GraphImplementation.EPS)
 
686
Show(convergenceGraphIS)