Hello everybody,
With one of my student we are in the process of building a PCE for a function of continuous and discrete distribution of parameters.
To put it simple we are trying to fit functions of the type (see provided script):
with
- x_1, x_2 uniform continuous variables
- x_3 and x_4 discrete variables taking the values -1 or +1 (each of them)
Our actual functions are more complex than this one.
One first build the base of the PCE then use an hyperbolic anisotropic enumeration function. The idea was to limit the order of the polynomials acting for x3 and x4 (per the documentation they are limited to order 1, Orthogonal polynomials — OpenTURNS 1.20 documentation, Krawtchouk polynomials from the Bernouilli distribution interpreted as a simple case of a Binomial distribution).
As one wants to catch the x_2^2 term we put a weight for the enumeration function. After some checks we have seen that the (infinite) basis along the x3, x4 axis have orthogonal polynoms up to order 1. Higher orders are present but with inf and nan values as expected from the recursion formula.
We modified the basis by forcing the enumeration to put in front the wanted terms with the getStrataCumulatedCardinal method up to a given order of the polynomials involving the discrete terms.
Then we applied the cleaning strategy. It is working but not everytime. Sometimes the solver indicates that we are trying to identifiy on an orthogonal basis with too high orders on the x3, x4 axis. We obtained this by playing with the number provided for the in the getStrataCumulatedCardinal method.
We tried, alternatively, to build the basis ourselves with chosen polynomials but it does not work as OT characterize our basis as non orthogonal. Even if we think that it is actual a proper basis.
Finally we build our metamodel as follow:
Building 4 functions for each combinations of (x3,x4) 2-plet: f(x1,x2,+1,+1)=f_++(x1,x2), etc …
Building 4 metamodels for (x1,x2) continuous variables
Regrouping everything, it works perfectly.
The solution is somewhat OK as it is a little bit inaesthetic, is it possible to perform the task without the trick of splitting our function in 4 pieces ?
Regards
import openturns as ot
import openturns.viewer as viewer
import numpy as np
inputs = ['x1', 'x2', 'x3', 'x4']
formula = ['x1 * (x2 - 1)^2 + x3 * x1 * x2 * x4']
model = ot.SymbolicFunction(inputs, formula)
distribution = ot.ComposedDistribution(
[ot.Uniform(1, 10), ot.Uniform(1, 10), 2*ot.Bernoulli(0.5)-1, 2*ot.Bernoulli(0.5)-1])
# Training sample
sample_size = 1000
X = ot.LHSExperiment(distribution, sample_size).generate()
Y = model(X)
# Sobol model
model_sa = ot.SaltelliSensitivityAlgorithm(distribution, sample_size, model, True)
# Multivariate basis
dim = distribution.getDimension()
uni_basis = ot.PolynomialFamilyCollection(dim)
for i in range(dim):
distrib = distribution.getMarginal(i)
uni_basis[i] = ot.StandardDistributionPolynomialFactory(distrib)
weight = 5
enumerate_function = ot.HyperbolicAnisotropicEnumerateFunction([1, 1, weight, weight], 1)
multi_basis = ot.OrthogonalProductPolynomialFactory(uni_basis, enumerate_function)
p = enumerate_function.getStrataCumulatedCardinal(
enumerate_function.inverse([1, 1, 1, 1]) + 1)
n = 9
truncature_strategy = ot.CleaningStrategy(multi_basis, p, n, 1e-3, True)
basis_factory = ot.LARS() # Least angle regression
fit_algo = ot.CorrectedLeaveOneOut()
approx_algo = ot.LeastSquaresMetaModelSelectionFactory(basis_factory, fit_algo)
coef_eval_strategy = ot.LeastSquaresStrategy(
ot.MonteCarloExperiment(sample_size), approx_algo)
# PCE metamodel
algo = ot.FunctionalChaosAlgorithm(
X, Y, distribution, truncature_strategy, coef_eval_strategy)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()
coefficients = result.getCoefficients()
for i, j in enumerate(result.getIndices()):
indices = enumerate_function(j)
degree = sum(indices)
print(str(int(degree)) + ' ' + str(indices) + ' ' + str(coefficients[i]))
# Sobol metamodel
chaosSI = ot.FunctionalChaosSobolIndices(result)
metamodel_sa = ot.SaltelliSensitivityAlgorithm(
distribution, sample_size, metamodel, True)
# Validation sample
X_val = distribution.getSample(sample_size)
Y_val = model(X_val)
# Validation
val = ot.MetaModelValidation(X_val, Y_val, metamodel)
Q2 = val.computePredictivityFactor()
graph = val.drawValidation()
graph.setTitle("Q² = %.2f%%" % (Q2*100))
view = viewer.View(graph)
graph = model_sa.draw()
graph.add(metamodel_sa.draw())
view = viewer.View(graph)
# fixed_index = [1,2,3]
# fixed_values = [2,1,0]
# model_parametric = ot.ParametricFunction(model, fixed_index, fixed_values)
# metamodel_parametric = ot.ParametricFunction(metamodel, fixed_index, fixed_values)
# graph = ot.Graph('', 'x1', 'y0', True, 'bottomright')
# model_curve = model_parametric.getMarginal(0).draw(0, 10).getDrawable(0)
# model_curve.setLegend('Model')
# metamodel_curve = metamodel_parametric.getMarginal(0).draw(0, 10).getDrawable(0)
# metamodel_curve.setLegend('Metamodel')
# metamodel_curve.setLineStyle('dashed')
# metamodel_curve.setColor('red')
# graph.add(model_curve)
# graph.add(metamodel_curve)
# viewer.View(graph)