Note

Go to the end to download the full example code.

Plot enumeration function

This example illustrates the enumeration functions which are ijections between \Nset into \Nset^\inputDim. Refer to Multivariate indices enumeration functions to get the precise description of the enumerate functions. We detail here the bijections:

  • Linear enumeration function

  • Hyperbolic enumeration function

  • Anisotropic hyperbolic enumeration function

  • Infinity norm enumeration function

These bijections are used in the in the functional chaos expansion setting. In order to build up a multivariate basis \{\psi_{\vect{\alpha}},\vect{\alpha} \in \Nset^\inputDim\} by tensorization of univariate basis terms, we need to enumerate the multi-indices \vect{\alpha} \in \Nset^\inputDim. In this example, we interprete the impact of the different enumeration functions within the functional chaos expansion setting.

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. If the basis is polynomial, then the length corresponds to the total degree of the polynomial. Within a strata, the multi-indices are ordered according to the “graded reverse lexicographic ordering” in [sullivan2015]. This is named the linear enumeration rule in the library.

We print the ordered elements in dimension 4 to illustrate the ordering of the multi-indices. We only consider the 4 first strata and we print their associated length (which are equal to 0 up to 3).

dim = 4
enum_func = ot.LinearEnumerateFunction(dim)
td_prev = 0
print("#  | multi-index     | length")
print("---+-----------------+-------------")
i_max = enum_func.getStrataCumulatedCardinal(3)
for i in range(i_max):
    multi_index = enum_func(i)
    td = sum(multi_index)
    if td > td_prev:
        td_prev = td
        print("---+-----------------+-------------")
    print(f"{i:2} |       {multi_index} |           {td}")

#  | multi-index     | length
---+-----------------+-------------
 0 |       [0,0,0,0] |           0
---+-----------------+-------------
 1 |       [1,0,0,0] |           1
 2 |       [0,1,0,0] |           1
 3 |       [0,0,1,0] |           1
 4 |       [0,0,0,1] |           1
---+-----------------+-------------
 5 |       [2,0,0,0] |           2
 6 |       [1,1,0,0] |           2
 7 |       [1,0,1,0] |           2
 8 |       [1,0,0,1] |           2
 9 |       [0,2,0,0] |           2
10 |       [0,1,1,0] |           2
11 |       [0,1,0,1] |           2
12 |       [0,0,2,0] |           2
13 |       [0,0,1,1] |           2
14 |       [0,0,0,2] |           2
---+-----------------+-------------
15 |       [3,0,0,0] |           3
16 |       [2,1,0,0] |           3
17 |       [2,0,1,0] |           3
18 |       [2,0,0,1] |           3
19 |       [1,2,0,0] |           3
20 |       [1,1,1,0] |           3
21 |       [1,1,0,1] |           3
22 |       [1,0,2,0] |           3
23 |       [1,0,1,1] |           3
24 |       [1,0,0,2] |           3
25 |       [0,3,0,0] |           3
26 |       [0,2,1,0] |           3
27 |       [0,2,0,1] |           3
28 |       [0,1,2,0] |           3
29 |       [0,1,1,1] |           3
30 |       [0,1,0,2] |           3
31 |       [0,0,3,0] |           3
32 |       [0,0,2,1] |           3
33 |       [0,0,1,2] |           3
34 |       [0,0,0,3] |           3

We define a function that plots the successive stratas with different colors.

def draw_stratas(enum_func, maximum_strata_index):
    """
    Plot enumeration rule by stratas.

    Parameters
    ----------
    enum_func : openturns.EnumerateFunction
        Enumerate function

    Returns
    -------
    graph : openturns.Graph
        Plot of the multi-indices colored by stratas
    """
    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()}")
    elif enum_func.__class__.__name__ == "NormInfEnumerateFunction":
        graph.setTitle("Infinity Norm enumeration rule")
    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

We consider the dimension 2 and we plot the first strata.

dim = 2
enum_func = ot.LinearEnumerateFunction(dim)
graph = draw_stratas(enum_func, 7)
view = otv.View(graph, axes_kw={"aspect": "equal"})
Linear enumeration rule

When the number of input dimensions of a functional chaos expansion (FCE) increases, each multi-index corresponds to a coefficient in the expansion. Hence, the number of multi-indices represents the number of coefficients in the FCE. We plot the number of terms in the basis depending on the maximum total length for several dimension values. We observe the exponential increase of the number of terms with the dimension \inputDim (curse of dimensionality).

graph = ot.Graph("Linear enumeration", "Total degree", "Number of coefficients")
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)})
Linear enumeration

The hyperbolic enumeration function is based on the q-norm. We plot the hyperbolic quasi norm for different values of q. With q = 1 (with isotropy), stratas are hyperplanes: we recover 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"})
Hyperbolic quasi norm, q = 1.0, q = 0.75, q = 0.5, q = 0.25

We plot the multi-indices of the hyperbolic isotropic enumeration rule by stratas. The lower the value of 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), 7)
)
grid.setGraph(
    0, 1, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.75), 7)
)
grid.setGraph(
    1, 0, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.5), 7)
)
grid.setGraph(
    1, 1, draw_stratas(ot.HyperbolicAnisotropicEnumerateFunction(dim, 0.25), 7)
)
grid.setTitle("Hyperbolic rule")
view = otv.View(grid, axes_kw={"aspect": "equal"})
Hyperbolic rule, q=1.0, q=0.75, q=0.5, q=0.25

Interaction multi-indices are presented in the center of the (\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 axes represent multivariate polynomials with high marginal degrees. When q is close to zero, these polynomials with high marginal degrees appear sooner with the hyperbolic enumeration rule.

We plot the number of terms in the basis depending on the maximum total degree in dimension d = 5 for several q -norm values. We observe that the gap between high and low values of 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)})
Hyperbolic enumeration. dim. = 5

We plot the multi-indices of the hyperbolic anisotropic enumeration rule by stratas. This enumerate function emphasizes multi-indices whose components are larger when the associated weights are smaller.

grid = ot.GridLayout(2, 2)
weights = [0.4, 0.6]
maximum_strata_index = 14
grid.setGraph(
    0,
    0,
    draw_stratas(
        ot.HyperbolicAnisotropicEnumerateFunction(weights, 1.0), maximum_strata_index
    ),
)
grid.setGraph(
    0,
    1,
    draw_stratas(
        ot.HyperbolicAnisotropicEnumerateFunction(weights, 0.7), maximum_strata_index
    ),
)
grid.setGraph(
    1,
    0,
    draw_stratas(
        ot.HyperbolicAnisotropicEnumerateFunction(weights, 0.5), maximum_strata_index
    ),
)
grid.setGraph(
    1,
    1,
    draw_stratas(
        ot.HyperbolicAnisotropicEnumerateFunction(weights, 0.25), maximum_strata_index
    ),
)
grid.setTitle("Hyperbolic anisotropic rule, weights = [0.4, 0.6]")
view = otv.View(grid, axes_kw={"aspect": "equal"})
Hyperbolic anisotropic rule, weights = [0.4, 0.6], q=1.0, q=0.7, q=0.5, q=0.25

Now we use the infinity norm enumeration function. We illustrate the enumeration in dimension 2. We plot the first stratas.

enum_func = ot.NormInfEnumerateFunction(2)
graph = draw_stratas(enum_func, 7)
view = otv.View(graph, axes_kw={"aspect": "equal"})
Infinity Norm enumeration rule

We print the 3 first stratas in dimension 3.

print("#  | multi-index ")
print("---+-------------")
enum_func = ot.NormInfEnumerateFunction(3)
i_max = enum_func.getStrataCumulatedCardinal(3)
for i in range(i_max):
    multi_index = enum_func(i)
    print(f"{i:2} |       {multi_index}")

#  | multi-index
---+-------------
 0 |       [0,0,0]
 1 |       [1,0,0]
 2 |       [0,1,0]
 3 |       [1,1,0]
 4 |       [0,0,1]
 5 |       [1,0,1]
 6 |       [0,1,1]
 7 |       [1,1,1]
 8 |       [2,0,0]
 9 |       [2,1,0]
10 |       [0,2,0]
11 |       [1,2,0]
12 |       [2,2,0]
13 |       [2,0,1]
14 |       [2,1,1]
15 |       [0,2,1]
16 |       [1,2,1]
17 |       [2,2,1]
18 |       [0,0,2]
19 |       [1,0,2]
20 |       [2,0,2]
21 |       [0,1,2]
22 |       [1,1,2]
23 |       [2,1,2]
24 |       [0,2,2]
25 |       [1,2,2]
26 |       [2,2,2]
27 |       [3,0,0]
28 |       [3,1,0]
29 |       [3,2,0]
30 |       [0,3,0]
31 |       [1,3,0]
32 |       [2,3,0]
33 |       [3,3,0]
34 |       [3,0,1]
35 |       [3,1,1]
36 |       [3,2,1]
37 |       [0,3,1]
38 |       [1,3,1]
39 |       [2,3,1]
40 |       [3,3,1]
41 |       [3,0,2]
42 |       [3,1,2]
43 |       [3,2,2]
44 |       [0,3,2]
45 |       [1,3,2]
46 |       [2,3,2]
47 |       [3,3,2]
48 |       [0,0,3]
49 |       [1,0,3]
50 |       [2,0,3]
51 |       [3,0,3]
52 |       [0,1,3]
53 |       [1,1,3]
54 |       [2,1,3]
55 |       [3,1,3]
56 |       [0,2,3]
57 |       [1,2,3]
58 |       [2,2,3]
59 |       [3,2,3]
60 |       [0,3,3]
61 |       [1,3,3]
62 |       [2,3,3]
63 |       [3,3,3]

Display the graphs.

otv.View.ShowAll()