Polynomial chaos exploitationΒΆ

In this example we are going to create a global approximation of a model response using functional chaos and expose the associated results:

  • the composed model: h: \underline{Z}^{\strut} \longrightarrow \underline{Y} = g \circ T^{-1}(\underline{Z}), which is the model of the reduced variables \underline{Z}. We have \displaystyle h = \sum_{k \in \mathbb N} \underline{\alpha}_k \Psi_k,

  • the coefficients of the polynomial approximation : (\underline{\alpha}_k)_{k \in K},

  • the composed meta model: \hat{h}, which is the model of the reduced variables reduced to the truncated multivariate basis (\Psi_k)_{k \in K}. We have \displaystyle \hat{h} = \sum_{k \in K} \underline{\alpha}_k \Psi_k,

  • the meta model: \displaystyle \hat{g} : \underline{X} \longrightarrow Y = \hat{h} \circ T(\underline{X}) which is the polynomial chaos approximation as a NumericalMathFunction. We have \displaystyle \hat{g} = \sum_{k \in K} \underline{\alpha}_k \Psi_k \circ T,

  • the truncated multivariate basis : (\Psi_k)_{k \in K},

  • the indices K,

  • the composition of each polynomial of the truncated multivariate basis \Psi_k,

  • the distribution \mu of the transformed variables \underline{Z},

[1]:
from __future__ import print_function
import openturns as ot
[2]:
# prepare some X/Y data
ot.RandomGenerator.SetSeed(0)
dimension = 2
input_names = ['x1', 'x2']
formulas = ['cos(x1 + x2)', '(x2 + 1) * exp(x1 - 2 * x2)']
model = ot.SymbolicFunction(input_names, formulas)
distribution = ot.Normal(dimension)
x  = distribution.getSample(30)
y = model(x)
[3]:
# create a functional chaos algorithm
algo = ot.FunctionalChaosAlgorithm(x, y)
algo.run()
[4]:
# Stream out the result
result = algo.getResult()
[5]:
# Get the polynomial chaos coefficients
# All the coefficients
result.getCoefficients()
[5]:
v0v1
00.425628570114641450.19814108237713063
10.00.470485758115012
2-0.099320894283325840.7885721936793412
3-0.246250852360200920.3232951290142421
4-0.19204843536289498-1.8462421285361612
50.01.4246957691984192
6-0.476436113471704840.0
70.0-1.029321781268946
80.0-0.1355602355790481
90.0-0.28170119639243424
10-0.02491833729348950.0
110.0-0.30064668963377356
120.075582424033229260.0
130.00.34014180081319967
140.188744950428186680.0
150.00.03466594712788605
[6]:
# The coefficients of marginal i
i = 1
result.getCoefficients()[i]
[6]:

[0,0.470486]

[7]:
# Get the indices of the selected polynomials : K
subsetK = result.getIndices()
subsetK
[7]:

[0,1,2,3,4,6,7,9,12,13,16,17,18,19,20,21]#16

[8]:
# Get the composition of the polynomials
# of the truncated multivariate basis
for i in range(subsetK.getSize()):
    print("Polynomial number ", i, " in truncated basis <-> polynomial number ",
          subsetK[i], " = ", ot.LinearEnumerateFunction(dimension)(subsetK[i]), " in complete basis")
Polynomial number  0  in truncated basis <-> polynomial number  0  =  [0,0]  in complete basis
Polynomial number  1  in truncated basis <-> polynomial number  1  =  [1,0]  in complete basis
Polynomial number  2  in truncated basis <-> polynomial number  2  =  [0,1]  in complete basis
Polynomial number  3  in truncated basis <-> polynomial number  3  =  [2,0]  in complete basis
Polynomial number  4  in truncated basis <-> polynomial number  4  =  [1,1]  in complete basis
Polynomial number  5  in truncated basis <-> polynomial number  6  =  [3,0]  in complete basis
Polynomial number  6  in truncated basis <-> polynomial number  7  =  [2,1]  in complete basis
Polynomial number  7  in truncated basis <-> polynomial number  9  =  [0,3]  in complete basis
Polynomial number  8  in truncated basis <-> polynomial number  12  =  [2,2]  in complete basis
Polynomial number  9  in truncated basis <-> polynomial number  13  =  [1,3]  in complete basis
Polynomial number  10  in truncated basis <-> polynomial number  16  =  [4,1]  in complete basis
Polynomial number  11  in truncated basis <-> polynomial number  17  =  [3,2]  in complete basis
Polynomial number  12  in truncated basis <-> polynomial number  18  =  [2,3]  in complete basis
Polynomial number  13  in truncated basis <-> polynomial number  19  =  [1,4]  in complete basis
Polynomial number  14  in truncated basis <-> polynomial number  20  =  [0,5]  in complete basis
Polynomial number  15  in truncated basis <-> polynomial number  21  =  [6,0]  in complete basis
[9]:
# Get the multivariate basis
# as a collection of Function
reduced = result.getReducedBasis()
[10]:
# Get the orthogonal basis
orthgBasis = result.getOrthogonalBasis()
orthgBasis
[10]:

class=OrthogonalProductPolynomialFactory univariate polynomial collection=[class=OrthogonalUniVariatePolynomialFamily implementation=class=StandardDistributionPolynomialFactory hasSpecificFamily=true specificFamily=class=OrthogonalUniVariatePolynomialFamily implementation=class=HistogramPolynomialFactory measure=class=Histogram name=Histogram dimension=1 first=-1 width=class=Point name=Unnamed dimension=5 values=[0.4,0.4,0.4,0.4,0.4] height=class=Point name=Unnamed dimension=5 values=[0.25,0.25,0.666667,1.08333,0.25] legendre=class=LegendreFactory measure=class=Uniform name=Uniform dimension=1 a=-1 b=1,class=OrthogonalUniVariatePolynomialFamily implementation=class=StandardDistributionPolynomialFactory hasSpecificFamily=false orthonormalization algorithm=class=OrthonormalizationAlgorithm implementation=class=AdaptiveStieltjesAlgorithm measure=class=Trapezoidal name=Trapezoidal dimension=1 a=-1 b=-0.530738 c=0.0480939 d=1 h=0.775545 monicRecurrenceCoefficients=[class=Point name=Unnamed dimension=3 values=[1,0.098496,0],class=Point name=Unnamed dimension=3 values=[1,0.0055215,-0.184989],class=Point name=Unnamed dimension=3 values=[1,-0.0074096,-0.222231],class=Point name=Unnamed dimension=3 values=[1,0.0105003,-0.237692],class=Point name=Unnamed dimension=3 values=[1,-0.00366255,-0.243879],class=Point name=Unnamed dimension=3 values=[1,0.000574247,-0.242444],class=Point name=Unnamed dimension=3 values=[1,0.000785537,-0.246691],class=Point name=Unnamed dimension=3 values=[1,0.000879275,-0.247369],class=Point name=Unnamed dimension=3 values=[1,-0.00173165,-0.247045]] monicSquaredNorms=class=Point name=Unnamed dimension=9 values=[1,0.184989,0.0411103,0.00977157,0.00238308,0.000577763,0.000142529,3.52572e-05,8.71013e-06] isElliptical=false] measure=class=ComposedDistribution name=ComposedDistribution dimension=2 copula=class=IndependentCopula name=IndependentCopula dimension=2 marginal[0]=class=Histogram name=Histogram dimension=1 first=-1 width=class=Point name=Unnamed dimension=5 values=[0.4,0.4,0.4,0.4,0.4] height=class=Point name=Unnamed dimension=5 values=[0.25,0.25,0.666667,1.08333,0.25] marginal[1]=class=Trapezoidal name=Trapezoidal dimension=1 a=-1 b=-0.530738 c=0.0480939 d=1 h=0.775545

[11]:
# Get the distribution of variables Z
orthgBasis.getMeasure()
[11]:

ComposedDistribution(Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25}), Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1), IndependentCopula(dimension = 2))

[12]:
# Get the composed model which is the model of the reduced variables Z
result.getComposedModel()
[12]:

(DatabaseEvaluation
input sample :
[ X0 X1 ]
0 : [ 0.608202 -1.26617 ]
1 : [ -0.438266 1.20548 ]
2 : [ -2.18139 0.350042 ]
3 : [ -0.355007 1.43725 ]
4 : [ 0.810668 0.793156 ]
5 : [ -0.470526 0.261018 ]
6 : [ -2.29006 -1.28289 ]
7 : [ -1.31178 -0.0907838 ]
8 : [ 0.995793 -0.139453 ]
9 : [ -0.560206 0.44549 ]
10 : [ 0.322925 0.445785 ]
11 : [ -1.03808 -0.856712 ]
12 : [ 0.473617 -0.125498 ]
13 : [ 0.351418 1.78236 ]
14 : [ 0.0702074 -0.781366 ]
15 : [ -0.721533 -0.241223 ]
16 : [ -1.78796 0.40136 ]
17 : [ 1.36783 1.00434 ]
18 : [ 0.741548 -0.0436123 ]
19 : [ 0.539345 0.29995 ]
20 : [ 0.407717 -0.485112 ]
21 : [ -0.382992 -0.752817 ]
22 : [ 0.257926 1.96876 ]
23 : [ -0.671291 1.85579 ]
24 : [ 0.0521593 0.790446 ]
25 : [ 0.716353 -0.743622 ]
26 : [ 0.184356 -1.53073 ]
27 : [ 0.655027 0.538071 ]
28 : [ 1.73821 -0.958722 ]
29 : [ 0.377922 -0.181004 ]
output sample :
[ y0 y1 ]
0 : [ 0.791234 -6.153 ]
1 : [ 0.719848 0.127674 ]
2 : [ -0.257609 0.075673 ]
3 : [ 0.46935 0.0964592 ]
4 : [ -0.0330217 0.825582 ]
5 : [ 0.978133 0.467366 ]
6 : [ -0.9084 -0.372691 ]
7 : [ 0.167439 0.293644 ]
8 : [ 0.655206 3.07871 ]
9 : [ 0.993427 0.338667 ]
10 : [ 0.718808 0.818737 ]
11 : [ -0.318354 0.28152 ]
12 : [ 0.940016 1.80491 ]
13 : [ -0.533709 0.111917 ]
14 : [ 0.757606 1.11916 ]
15 : [ 0.571259 0.59742 ]
16 : [ 0.183152 0.105058 ]
17 : [ -0.718312 1.05597 ]
18 : [ 0.76617 2.19061 ]
19 : [ 0.667988 1.22357 ]
20 : [ 0.997007 2.04242 ]
21 : [ 0.421399 0.759585 ]
22 : [ -0.609865 0.0749114 ]
23 : [ 0.376759 0.0356671 ]
24 : [ 0.665521 0.388187 ]
25 : [ 0.999628 2.32215 ]
26 : [ 0.222539 -13.6308 ]
27 : [ 0.368781 1.00946 ]
28 : [ 0.711272 1.59716 ]
29 : [ 0.980674 1.71644 ])o(((| y0 = Normal(mu = 0, sigma = 1) -> y0 : Histogram(origin = -2.29006, {w0 = 0.805655, h0 = 0.124123}, {w1 = 0.805655, h1 = 0.124123}, {w2 = 0.805655, h2 = 0.330994}, {w3 = 0.805655, h3 = 0.537865}, {w4 = 0.805655, h4 = 0.124123})
| y1 = Normal(mu = 0, sigma = 1) -> y1 : Trapezoidal(a = -1.75839, b = -0.781359, c = 0.423804, d = 2.40573)
)o(class=LinearEvaluation name=Unnamed center=[0,0] constant=[0,0] linear=[[ 1 0.0108329 ]
[ 0 0.999941 ]]))o(| y0 = Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25}) -> y0 : Normal(mu = 0, sigma = 1)
| y1 = Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1) -> y1 : Normal(mu = 0, sigma = 1)
))

[13]:
# Get the composed meta model which is the model of the reduced variables Z
# within the reduced polynomials basis
result.getComposedMetaModel()
[13]:

[0.425629,0.198141] + [0,0.470486] * (-0.29173 + 2.18797 * x0) + [-0.0993209,0.788572] * (0.229006 + 2.32502 * x1) + [-0.246251,0.323295] * (-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) + [-0.192048,-1.84624] * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2) + [0,1.4247] * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3) + [-0.476436,0] * ((-0.29173 + 2.18797 * x0) * (0.229006 + 2.32502 * x1)) + [0,-1.02932] * (0.893658 - 1.0145 * x1 - 13.1947 * x1^2 + 2.19409 * x1^3 + 20.4848 * x1^4) + [0,-0.13556] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (0.229006 + 2.32502 * x1)) + [0,-0.281701] * ((-0.29173 + 2.18797 * x0) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [-0.0249183,0] * (0.284385 - 7.63723 * x0 - 4.67858 * x0^2 + 78.3356 * x0^3 + 9.59256 * x0^4 - 185.818 * x0^5 - 5.17002 * x0^6 + 120.521 * x0^7) + [0,-0.300647] * (-0.202696 - 7.72688 * x1 + 5.82995 * x1^2 + 87.3021 * x1^3 - 20.8178 * x1^4 - 231.97 * x1^5 + 17.6506 * x1^6 + 168.413 * x1^7) + [0.0755824,0] * ((0.400359 - 4.06203 * x0 - 0.288198 * x0^2 + 7.74638 * x0^3) * (0.229006 + 2.32502 * x1)) + [0,0.340142] * ((-0.29173 + 2.18797 * x0) * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3)) + [0.188745,0] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [0,0.0346659] * (0.95105 + 2.92377 * x0 - 42.0158 * x0^2 - 17.3514 * x0^3 + 250.399 * x0^4 + 29.114 * x0^5 - 454.203 * x0^6 - 14.8305 * x0^7 + 250.574 * x0^8)

[14]:
# Get the meta model which is the composed meta model combined with the
# iso probabilistic transformation
result.getMetaModel()
[14]:

([0.425629,0.198141] + [0,0.470486] * (-0.29173 + 2.18797 * x0) + [-0.0993209,0.788572] * (0.229006 + 2.32502 * x1) + [-0.246251,0.323295] * (-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) + [-0.192048,-1.84624] * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2) + [0,1.4247] * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3) + [-0.476436,0] * ((-0.29173 + 2.18797 * x0) * (0.229006 + 2.32502 * x1)) + [0,-1.02932] * (0.893658 - 1.0145 * x1 - 13.1947 * x1^2 + 2.19409 * x1^3 + 20.4848 * x1^4) + [0,-0.13556] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (0.229006 + 2.32502 * x1)) + [0,-0.281701] * ((-0.29173 + 2.18797 * x0) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [-0.0249183,0] * (0.284385 - 7.63723 * x0 - 4.67858 * x0^2 + 78.3356 * x0^3 + 9.59256 * x0^4 - 185.818 * x0^5 - 5.17002 * x0^6 + 120.521 * x0^7) + [0,-0.300647] * (-0.202696 - 7.72688 * x1 + 5.82995 * x1^2 + 87.3021 * x1^3 - 20.8178 * x1^4 - 231.97 * x1^5 + 17.6506 * x1^6 + 168.413 * x1^7) + [0.0755824,0] * ((0.400359 - 4.06203 * x0 - 0.288198 * x0^2 + 7.74638 * x0^3) * (0.229006 + 2.32502 * x1)) + [0,0.340142] * ((-0.29173 + 2.18797 * x0) * (-0.207607 - 4.12182 * x1 + 0.977306 * x1^2 + 10.1162 * x1^3)) + [0.188745,0] * ((-0.922997 + 0.0686243 * x0 + 4.03168 * x0^2) * (-0.909686 + 0.513016 * x1 + 4.93202 * x1^2)) + [0,0.0346659] * (0.95105 + 2.92377 * x0 - 42.0158 * x0^2 - 17.3514 * x0^3 + 250.399 * x0^4 + 29.114 * x0^5 - 454.203 * x0^6 - 14.8305 * x0^7 + 250.574 * x0^8))o((| y0 = Normal(mu = 0, sigma = 1) -> y0 : Histogram(origin = -1, {w0 = 0.4, h0 = 0.25}, {w1 = 0.4, h1 = 0.25}, {w2 = 0.4, h2 = 0.666667}, {w3 = 0.4, h3 = 1.08333}, {w4 = 0.4, h4 = 0.25})
| y1 = Normal(mu = 0, sigma = 1) -> y1 : Trapezoidal(a = -1, b = -0.530738, c = 0.0480939, d = 1)
)o((class=LinearEvaluation name=Unnamed center=[0,0] constant=[0,0] linear=[[ 1 -0.0108335 ]
[ 0 1.00006 ]])o(| y0 = Histogram(origin = -2.29006, {w0 = 0.805655, h0 = 0.124123}, {w1 = 0.805655, h1 = 0.124123}, {w2 = 0.805655, h2 = 0.330994}, {w3 = 0.805655, h3 = 0.537865}, {w4 = 0.805655, h4 = 0.124123}) -> y0 : Normal(mu = 0, sigma = 1)
| y1 = Trapezoidal(a = -1.75839, b = -0.781359, c = 0.423804, d = 2.40573) -> y1 : Normal(mu = 0, sigma = 1)
)))

[15]:
# Get the projection strategy
algo.getProjectionStrategy()
[15]:

class=ProjectionStrategy implementation=class=LeastSquaresStrategy experiment=class=FixedExperiment name=Unnamed sample=class=Sample name=Normal implementation=class=SampleImplementation name=Normal size=30 dimension=2 description=[X0,X1] data=[[0.608202,-1.26617],[-0.438266,1.20548],[-2.18139,0.350042],[-0.355007,1.43725],[0.810668,0.793156],[-0.470526,0.261018],[-2.29006,-1.28289],[-1.31178,-0.0907838],[0.995793,-0.139453],[-0.560206,0.44549],[0.322925,0.445785],[-1.03808,-0.856712],[0.473617,-0.125498],[0.351418,1.78236],[0.0702074,-0.781366],[-0.721533,-0.241223],[-1.78796,0.40136],[1.36783,1.00434],[0.741548,-0.0436123],[0.539345,0.29995],[0.407717,-0.485112],[-0.382992,-0.752817],[0.257926,1.96876],[-0.671291,1.85579],[0.0521593,0.790446],[0.716353,-0.743622],[0.184356,-1.53073],[0.655027,0.538071],[1.73821,-0.958722],[0.377922,-0.181004]] weights=class=Point name=Unnamed dimension=30 values=[0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333,0.0333333]