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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
"""
Generate low discrepancy sequences
==================================
"""
# %%
# In this examples we are going to expose the available low discrepancy sequences in order to approximate some integrals.
#
# The following low-discrepancy sequences are available:
#
# - Sobol
# - Faure
# - Halton
# - reverse Halton
# - Haselgrove
#
# To illustrate these sequences we generate their first 1024 points and compare with the sequence obtained from
# the pseudo random generator (Merse Twister) as the latter has a higher discrepancy.
# %%
import openturns as ot
import openturns.viewer as otv
# %%
# 1. Sobol sequence
dimension = 2
size = 1024
sequence = ot.SobolSequence(dimension)
sample = sequence.generate(size)
graph = ot.Graph("", "", "", True, "")
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# 2. Halton sequence
dimension = 2
sequence = ot.HaltonSequence(dimension)
sample = sequence.generate(size)
graph = ot.Graph("Halton", "", "", True, "")
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# 3. Halton sequence in high dimension: bad filling in upper dimensions
dimension = 20
sequence = ot.HaltonSequence(dimension)
sample = sequence.generate(size).getMarginal([dimension - 2, dimension - 1])
graph = ot.Graph(
"Halton (" + str(dimension - 2) + "," + str(dimension - 1) + ")",
"dim " + str(dimension - 2),
"dim " + str(dimension - 1),
True,
"",
)
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# 4. Scrambled Halton sequence in high dimension
dimension = 20
sequence = ot.HaltonSequence(dimension)
sequence.setScrambling("RANDOM")
sample = sequence.generate(size).getMarginal([dimension - 2, dimension - 1])
graph = ot.Graph(
"Halton (" + str(dimension - 2) + "," + str(dimension - 1) + ")",
"dim " + str(dimension - 2),
"dim " + str(dimension - 1),
True,
"",
)
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# 5. Reverse Halton sequence
dimension = 2
sequence = ot.ReverseHaltonSequence(dimension)
sample = sequence.generate(size)
print(
"discrepancy=",
ot.LowDiscrepancySequenceImplementation.ComputeStarDiscrepancy(sample),
)
graph = ot.Graph("Reverse Halton", "", "", True, "")
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# 6. Haselgrove sequence
dimension = 2
sequence = ot.HaselgroveSequence(dimension)
sample = sequence.generate(size)
graph = ot.Graph("Haselgrove", "", "", True, "")
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# Compare with uniform random sequence
distribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * 2)
sample = distribution.getSample(size)
print(
"discrepancy=",
ot.LowDiscrepancySequenceImplementation.ComputeStarDiscrepancy(sample),
)
graph = ot.Graph("Mersenne Twister", "", "", True, "")
cloud = ot.Cloud(sample)
graph.add(cloud)
view = otv.View(graph)
# %%
# Display all figures
otv.View.ShowAll()
|
