1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
#! /usr/bin/env python
import sys
import openturns as ot
ot.TESTPREAMBLE()
ot.Log.Show(ot.Log.INFO)
def progress(percent):
sys.stderr.write("-- progress=" + str(percent) + "%\n")
def stop():
sys.stderr.write("-- stop?\n")
return False
dimension = 3
# gsobol
# meanTh = 1.0
# a = [0.0] * dimension
# inputVariables = [''] * dimension
# formula = '1.0'
# covTh = 1.0
# for i in range(dimension):
# a[i] = 0.5 * i
# covTh = covTh * (1.0 + 1.0 / (3.0 * (1.0 + a[i]) ** 2))
# inputVariables[i] = "xi" + str(i)
# formula += " * ((abs(4.0 * xi" + str(i) + " - 2.0) + " + str(a[i]) + ") / (1.0 + " + str(a[i]) + "))"
# covTh = covTh - 1.0
# model = ot.SymbolicFunction(inputVariables, [formula])
# distribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
# ishigami
formula = "sin(pi_*X1)+7*sin(pi_*X2)*sin(pi_*X2)+0.1*((pi_*X3)*(pi_*X3)*(pi_*X3)*(pi_*X3))*sin(pi_*X1)"
model = ot.SymbolicFunction(["X1", "X2", "X3"], [formula])
distribution = ot.JointDistribution([ot.Uniform(-1.0, 1.0)] * dimension)
estimator = ot.SaltelliSensitivityAlgorithm()
estimator.setUseAsymptoticDistribution(True)
algo = ot.SobolSimulationAlgorithm(distribution, model, estimator)
algo.setMaximumOuterSampling(250)
# size N of Sobol experiment at each iteration, total size is N*(d+2)
algo.setExperimentSize(10000)
algo.setBlockSize(97)
algo.setIndexQuantileLevel(0.05)
algo.setIndexQuantileEpsilon(1e-2)
# algo.setProgressCallback(progress)
# algo.setStopCallback(stop)
print("algo=", algo)
# Perform the simulation
algo.run()
# Stream out the result
result = algo.getResult()
print("result=", result)
print("FO=", result.getFirstOrderIndicesEstimate())
print("TO=", result.getTotalOrderIndicesEstimate())
foDist = result.getFirstOrderIndicesDistribution()
print(foDist)
toDist = result.getTotalOrderIndicesDistribution()
print(toDist)
convergenceGraph = algo.drawFirstOrderIndexConvergence()
convergenceGraph = algo.drawTotalOrderIndexConvergence()
graph = result.draw()
print(graph)
# from openturns.viewer import View
# View(convergenceGraph).ShowAll()
|
