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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
"""
Plot enumeration rules
----------------------
"""
# %%
# In order to build up a functional chaos multivariate basis :math:`\{\psi_{\idx},\idx \in \NM\}`
# by tensorization of univariate basis terms, we need to enumerate the multi-indices :math:`\vect{\alpha} \in \mathbb{N}^{n_X}`.
# In this example we are going to explore properties of these enumeration rules.
# Refer also to :ref:`enumeration_strategy` in the theoric documentation.
import openturns as ot
import openturns.viewer as otv
import math as m
# %%
# The simplest way to generate the multi-indices is to enumerate the terms of increasing length.
# In other words, we enumerate the multi-indices with length equal to 0, then 1, 2, 3, etc.
# This is called "graded reverse-lexicographic ordering" in [sullivan2015]_.
# This is named the linear enumeration rule in the library; let us instantiate it in the 2-dimensional case.
dim = 2
enum_func = ot.LinearEnumerateFunction(dim)
# %%
# Print the 25 first multi-indices and their associated length.
td_prev = 0
print("# | multi-index | length")
print("---+-------------+-------------")
for i in range(25):
multi_index = enum_func(i)
td = sum(multi_index)
if td > td_prev:
td_prev = td
print("---+-------------+-------------")
print(f"{i:2} | {multi_index} | {td}")
# %%
# Plot the multi-indices of the enumeration rule by stratas.
# In the specific case of the linear enumerate function, each strata contains
# multi-indices of identical length that is also the index of the strata.
def draw_stratas(enum_func):
"""
Plot enumeration rule by stratas.
Parameters
----------
enum_func : openturns.EnumerateFunction
Enumerate function
Returns
-------
graph : openturns.Graph
Plot of the multi-indices colored by stratas
"""
maximum_strata_index = 7
graph = ot.Graph("", "$\\alpha_1$", "$\\alpha_2$", True)
if enum_func.__class__.__name__ == "LinearEnumerateFunction":
graph.setTitle("Linear enumeration rule")
elif enum_func.__class__.__name__ == "HyperbolicAnisotropicEnumerateFunction":
graph.setTitle(f"q={enum_func.getQ()}")
offset = 0
for strata_index in range(maximum_strata_index):
strata_cardinal = enum_func.getStrataCardinal(strata_index)
sample_in_layer = [enum_func(idx + offset) for idx in range(strata_cardinal)]
offset += strata_cardinal
cloud = ot.Cloud(sample_in_layer)
cloud.setLegend(str(strata_index))
cloud.setPointStyle("circle")
graph.add(cloud)
graph.setIntegerXTick(True)
graph.setIntegerYTick(True)
graph.setLegendPosition("upper right")
return graph
graph = draw_stratas(enum_func)
view = otv.View(graph, axes_kw={"aspect": "equal"})
# %%
# When the number of input dimensions of a polynomial chaos expansion (PCE) increases,
# each multi-index corresponds to a coefficient in the expansion.
# Hence, the number of multi-indices represents the number of coefficients in the PCE.
# Plot the number of terms in the basis depending on the maximum total degree
# for several dimension values.
# We observe the exponential increase of the number of terms with the dimension
# :math:`d` (curse of dimensionality).
graph = ot.Graph("Linear enumeration", "Total degree", "Number of coefficients", True)
degree_maximum = 10
list_of_dimensions = [1, 5, 10, 15, 20]
point_styles = ["bullet", "circle", "fdiamond", "fsquare", "triangleup"]
for i in range(len(list_of_dimensions)):
dimension = list_of_dimensions[i]
number_of_coeff_array = ot.Sample(degree_maximum, 2)
for degree in range(1, 1 + degree_maximum):
number_of_coeff_array[degree - 1, 0] = degree
number_of_coeff_array[degree - 1, 1] = m.comb(degree + dimension, degree)
cloud = ot.Cloud(number_of_coeff_array)
cloud.setPointStyle(point_styles[i])
cloud.setLegend(f"dim.={dimension}")
graph.add(cloud)
graph.setLegendPosition("upper left")
graph.setIntegerXTick(True)
graph.setLogScale(ot.GraphImplementation.LOGY)
view = otv.View(graph, figure_kw={"figsize": (5, 4)})
# %%
# Plot the hyperbolic quasi norm for different values of :math:`q`.
# With :math:`q = 1` stratas are hyperplanes, and in case of isotropy
# it is equivalent to the linear enumeration rule.
def draw_qnorm(q):
def qnorm(x):
norm = 0.0
for xi in x:
norm += xi**q
norm = norm ** (1.0 / q)
return [norm]
f = ot.PythonFunction(2, 1, qnorm)
f.setInputDescription(["x1", "x2"])
graph = f.draw([0.0] * 2, [1.0] * 2)
graph.setTitle(f"q = {q}")
return graph
dln = ot.ResourceMap.GetAsUnsignedInteger("Contour-DefaultLevelsNumber")
ot.ResourceMap.SetAsUnsignedInteger("Contour-DefaultLevelsNumber", 5)
grid = ot.GridLayout(2, 2)
grid.setGraph(0, 0, draw_qnorm(1.0))
grid.setGraph(0, 1, draw_qnorm(0.75))
grid.setGraph(1, 0, draw_qnorm(0.5))
grid.setGraph(1, 1, draw_qnorm(0.25))
ot.ResourceMap.SetAsUnsignedInteger("Contour-DefaultLevelsNumber", dln)
grid.setTitle("Hyperbolic quasi norm")
view = otv.View(grid, axes_kw={"aspect": "equal"})
# %%
# Plot the multi-indices of the linear enumeration rule by stratas.
# The lower the value of :math:`q` the lower the number of interactions terms in stratas.
grid = ot.GridLayout(2, 2)
grid.setGraph(0, 0, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 1.0)))
grid.setGraph(0, 1, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.75)))
grid.setGraph(1, 0, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.5)))
grid.setGraph(1, 1, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.25)))
grid.setTitle("Hyperbolic rule")
view = otv.View(grid, axes_kw={"aspect": "equal"})
# %%
# Interaction multi-indices are presented in the center of the :math:`(\alpha_1, \alpha_2)` space.
# We see that when the quasi-norm parameter is close to zero, the enumeration rule
# shows less interaction multi-indices.
# Instead, multi-indices close to the :math:`x` and :math:`y` axes represent multivariate polynomials
# with high marginal degrees.
# When :math:`q` is close to zero, these polynomials with high marginal degrees appear
# sooner with the hyperbolic enumeration rule.
# %%
# Plot the number of terms in the basis depending on the maximum total degree
# in dimension :math:`d = 5` for several :math:`q` -norm values.
# We observe that the gap between high and low values of :math:`q` can lead to a gap
# in the numbers of terms of an order of magnitude.
dim = 5
graph = ot.Graph(
"Hyperbolic enumeration. dim. = %d" % (dim),
"Total degree",
"Number of coefficients",
True,
)
degree_maximum = 10
q_list = [0.2, 0.4, 0.6, 0.8, 1.0]
point_styles = ["bullet", "circle", "fdiamond", "fsquare", "triangleup"]
for i in range(len(q_list)):
q = q_list[i]
enum_func = ot.HyperbolicAnisotropicEnumerateFunction(dim, q)
number_of_coeff_array = ot.Sample(degree_maximum, 2)
for p in range(1, 1 + degree_maximum):
number_of_coeff_array[p - 1, 0] = p
number_of_coeff_array[p - 1, 1] = enum_func.getMaximumDegreeCardinal(p)
cloud = ot.Cloud(number_of_coeff_array)
cloud.setPointStyle(point_styles[i])
cloud.setLegend(f"$q={q}$")
graph.add(cloud)
graph.setLegendPosition("upper left")
graph.setIntegerXTick(True)
graph.setLogScale(ot.GraphImplementation.LOGY)
view = otv.View(graph, figure_kw={"figsize": (5, 4)})
otv.View.ShowAll()
# %%
# When the quasi-norm parameter is close to 1, then the hyperbolic rule is equal to the
# linear enumeration rule and the number of coefficients is larger.
#
# In practice, we often test several values of the parameter :math:`q`, in the :math:`[0.5, 0.9]` range,
# for example :math:`q \in \{0.5, 0.6, 0.7, 0.8, 0.9\}`.
# %%
# Reset default settings
ot.ResourceMap.Reload()
|
