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
|
#! /usr/bin/env python
import openturns as ot
import openturns.testing as ott
import openturns.experimental as otexp
ot.TESTPREAMBLE()
# Check the parameter constructor
mean = [0.0] * 2
std = [1.0] * 2
R = ot.CorrelationMatrix([[1.0, 0.8], [0.8, 1.0]])
distribution = ot.Normal(mean, std, R)
# These functions are obtained by using the Gram-Schmidt algorithm to make orthonormal
# wrt the distribution the first 2D tensorized Hermite polynomials
f0 = ot.SymbolicFunction(["x0", "x1"], ["1.0"])
f1 = ot.SymbolicFunction(["x0", "x1"], ["x0"])
f2 = ot.SymbolicFunction(["x0", "x1"], ["-4.0 * x0 / 3.0 + 5.0 * x1 / 3.0"])
f3 = ot.SymbolicFunction(["x0", "x1"], ["-1.0 / sqrt(2.0) + x0^2 / sqrt(2)"])
f4 = ot.SymbolicFunction(
["x0", "x1"],
[
"-4.0 / 3.0 - 1.885618083165693 * (-1.0 / sqrt(2.0) + x0^2 / sqrt(2.0)) + 5.0 * x0 * x1 / 3.0"
],
)
f5 = ot.SymbolicFunction(
["x0", "x1"],
[
"2.5141574442076222 + 16.0 / 9.0 * (-1.0 / sqrt(2.0) + x0^2 / sqrt(2.0)) - 3.142696805266918 * x0 * x1 + 25.0 / 9.0 * (-1.0 / sqrt(2.0) + x1^2 / sqrt(2.0))"
],
)
initialBasis = [f0, f1, f2, f3, f4, f5]
factory = otexp.FiniteOrthogonalFunctionFactory(initialBasis, distribution)
x = [0.5] * 2
kMax = len(initialBasis)
distribution = factory.getMeasure()
initialBasis = factory.getFunctionsCollection()
factory.setFunctionsCollection(initialBasis)
print("Factory=", factory)
print("Factory=", repr(factory))
functions = list()
for k in range(kMax):
fk = factory.build(k)
functions.append(fk)
value = fk(x)
print(f"FiniteOrthogonalFunction_{k}({x})={value}")
# Check orthonormality
M = ot.SymmetricMatrix(kMax)
integrationAlgo = ot.GaussLegendre([48] * 2)
for m in range(kMax):
for n in range(m + 1):
def wrapper(x):
return functions[m](x) * functions[n](x)[0] * distribution.computePDF(x)
kernel = ot.PythonFunction(distribution.getDimension(), 1, wrapper)
value = integrationAlgo.integrate(kernel, distribution.getRange())[0]
if abs(value) >= 1.0e-6:
M[m, n] = value
ott.assert_almost_equal(M, ot.IdentityMatrix(kMax))
|
