Advanced polynomial chaos constructionΒΆ

In this example we are going to expose advanced elements in the construction of a polynomial chaos algorithm:

  • Construction of the multivariate orthonormal basis

  • Truncature strategy of the multivariate orthonormal basis

  • Evaluation strategy of the approximation coefficients

[13]:
from __future__ import print_function
import openturns as ot
[14]:

# Create a function R^n --> R^p
# For example R^4 --> R
myModel = ot.SymbolicFunction(['x1', 'x2', 'x3', 'x4'], ['1+x1*x2 + 2*x3^2+x4^4'])

# Create a distribution of dimension n
# for example n=3 with indpendent components
Xdist = ot.ComposedDistribution([ot.Normal(), ot.Uniform(), ot.Gamma(2.75, 1.0), ot.Beta(2.5, 1.0, -1.0, 2.0)])

# Dimension of the input random vector
dim = 4
[15]:
#############################################################
# STEP 1: Construction of the multivariate orthonormal basis
#############################################################

# Create the univariate polynomial family collection
# which regroups the polynomial families for each direction
polyColl = ot.PolynomialFamilyCollection(dim)
# For information, with the Krawtchouk and Charlier families :
polyColl[0] = ot.KrawtchoukFactory()
polyColl[1] = ot.CharlierFactory()
# for information, with the automatic selection
for i in range(Xdist.getDimension()):
    polyColl[i] = ot.StandardDistributionPolynomialFactory(Xdist.getMarginal(i))
# We overload these factories as they are dedicated to discrete distributions
polyColl[0] = ot.HermiteFactory()
polyColl[1] = ot.LegendreFactory()
polyColl[2] = ot.LaguerreFactory(2.75)
# Parameter for the Jacobi factory : 'Probabilty' encoded with 1
polyColl[3] = ot.JacobiFactory(2.5, 3.5, 1)
[16]:
# Create the enumeration function
# LinearEnumerateFunction
enumerateFunction = ot.LinearEnumerateFunction(dim)

# HyperbolicAnisotropicEnumerateFunction
q = 0.4
enumerateFunction_1 = ot.HyperbolicAnisotropicEnumerateFunction(dim, q)
[17]:
# Create the multivariate orthonormal basis
# which is the the cartesian product of the univariate basis
multivariateBasis = ot.OrthogonalProductPolynomialFactory(
    polyColl, enumerateFunction)
[18]:
# Ask how many polynomials have degrees equal to k=5
k = 5
enumerateFunction.getStrataCardinal(k)
[18]:
56
[19]:
# Ask how many polynomials have degrees inferior or equal to k=5
enumerateFunction.getStrataCumulatedCardinal(k)
[19]:
126
[20]:
# Give the k-th term of the multivariate basis
# To calculate its degree, add the integers
k = 5
enumerateFunction(k)
[20]:

[2,0,0,0]

[21]:
# Build a term of the basis as a NumericalMathFunction
# Generally, we do not need to construct manually any term,
# all terms are constructed automatically by a strategy of construction of
# the basis
i = 5
Psi_i = multivariateBasis.build(i)
Psi_i
[21]:

-0.707107 + 0.707107 * x0^2

[22]:
# Get the measure mu associated to the multivariate basis
distributionMu = multivariateBasis.getMeasure()
distributionMu
[22]:

ComposedDistribution(Normal(mu = 0, sigma = 1), Uniform(a = -1, b = 1), Gamma(k = 3.75, lambda = 1, gamma = 0), Beta(r = 2.5, t = 3.5, a = -1, b = 1), IndependentCopula(dimension = 4))

[23]:
####################################################################
# STEP 2: Truncature strategy of the multivariate orthonormal basis
#############################################################
# FixedStrategy :
# all the polynomials af degree <=2
# which corresponds to the 15 first ones
p = 15
truncatureBasisStrategy = ot.FixedStrategy(multivariateBasis, p)

# SequentialStrategy :
# among the maximumCardinalBasis = 100 first polynomials
# of the multivariate basis those verfying the convergence criterion,
maximumCardinalBasis = 100
truncatureBasisStrategy_1 = ot.SequentialStrategy(
    multivariateBasis, maximumCardinalBasis)

# CleaningStrategy :
# among the maximumConsideredTerms = 500 first polynomials,
# those which have the mostSignificant = 50 most significant contributions
# with significance criterion significanceFactor = 10^(-4)
# The True boolean indicates if we are interested
# in the online monitoring of the current basis update
# (removed or added coefficients)
maximumConsideredTerms = 500
mostSignificant = 50
significanceFactor = 1.0e-4
truncatureBasisStrategy_2 = ot.CleaningStrategy(
    multivariateBasis, maximumConsideredTerms, mostSignificant, significanceFactor, True)
[24]:
################################################################
# STEP 3: Evaluation strategy of the approximation coefficients
#############################################################
# The technique illustrated is the Least Squares technique
# where the points come from an design of experiments
# Here : the Monte Carlo sampling technique
sampleSize = 100
evaluationCoeffStrategy = ot.LeastSquaresStrategy(
    ot.MonteCarloExperiment(sampleSize))

# You can specify the approximation algorithm
# This is the algorithm that generates a sequence of basis using Least
# Angle Regression
basisSequenceFactory = ot.LARS()

# This algorithm estimates the empirical error on each sub-basis using
# Leave One Out strategy
fittingAlgorithm = ot.CorrectedLeaveOneOut()
# Finally the metamodel selection algorithm embbeded in LeastSquaresStrategy
approximationAlgorithm = ot.LeastSquaresMetaModelSelectionFactory(
    basisSequenceFactory, fittingAlgorithm)
evaluationCoeffStrategy_2 = ot.LeastSquaresStrategy(
    ot.MonteCarloExperiment(sampleSize), approximationAlgorithm)

# Try integration
marginalDegrees = [2] * dim
evaluationCoeffStrategy_3 = ot.IntegrationStrategy(
    ot.GaussProductExperiment(distributionMu, marginalDegrees))
[25]:
#####################################################
# STEP 4: Creation of the Functional Chaos Algorithm
#############################################################
# FunctionalChaosAlgorithm :
# combination of the model : myModel
# the distribution of the input random vector : Xdist
# the truncature strategy of the multivariate basis
# and the evaluation strategy of the coefficients
polynomialChaosAlgorithm = ot.FunctionalChaosAlgorithm(
    myModel, Xdist, truncatureBasisStrategy, evaluationCoeffStrategy)