Question about metamodels choice

Methodology
1

Hello

I am currently building with openturns chaos polynomial and kriging metamodels
I have 20 input parameters and 3 outputs (scalars)

I build metamodel for each output independantly
with a sample of 200 computations (LHS)
I compute predictivity factor Q2 based on a test sample of 20 points
I could do more simulations if needed

for the first output it’s very good, for both type of metamodels (Q2 of 0.97 - 0.99)
but for the 2 other outputs it’s very bad (Q2 < 0.5, sometimes Q2 near 0 and < 0)

I have the feeling that these 2 results could be discontinuous (and physically it’s the case I believe)
one of them is zero in many cases (when geometric parameters are not too bad) and sometimes in some geometric configuration this result is not zero

what could I do in this context ? a better choice of some algorithm used by metamodel building process ?

Thanks in advance for your help !

Flore

2

Hi Flore,
In that case, you may test the following approach:

  • first classify the data,
  • then create a metamodel for each class.

This makes the task of the surrogate much simpler.
The overall metamodel is a composition :

Y = metamodel(classifier(X))

To perform classification, I have used logistic regression (and it improved the situation a lot), but other methods may perform better in some cases.
Best regards,
Michaël

3

Hello Flore,
In case you have not done it already, I would first suggest trying to visualize the distributions of the training output data (e.g., as histograms) for the 2 outputs which are poorly modeled. If the phenomena you are trying to model are indeed discontinuous, there is a good chance you’ll be able to see it right away.
Nevertheless, most of the standard surrogate modeling techniques present in the literature are unfortunately not particularly fit to properly model non-stationary functions. The most straightforward solution would be to follow Michaël’s suggestion and create separate surrogate models for different parts of the input design space, depending on the output value.

4

Hi Flore,

You may be interested by the otmixmod OpenTURNS module, in particular you should check this file. It implements a classification of the joint input/output space using a mixture on Gaussian distributions, then it build a polynomial chaos approximation for each classes and put everything into an ExpertMixture meta-model.

Don’t hesitate to contact me directly if you need further help!

Cheers

Régis

5

Thanks all!

6

histo_pdf_P54_sample_doe2
histo_pdf_P54_sample_doe21006×802 28.1 KB

histo_pdf_P55_sample_doe2
histo_pdf_P55_sample_doe21036×787 27.2 KB

histo_pdf_P57_sample_doe2
histo_pdf_P57_sample_doe21019×855 25.5 KB

I joined histograms
P54 is the result where Q2 is very good (near 1.0), it looks like a normal distrib or almost
P55 and P57 are the 2 other results
when P55 > 0.4 - 0.5 : this result have no more physical signification, even if a value is computed
and in many cases when P57 is big, corresponding R95 has no physical signification

7

I know nothing about classification, however I understand principles…
Define different metamodels from different input domains…
But I have 20 inputs…
Then ExpertMixture could be very intersting
I’ll take a look on that, thanks

8

I found some pdf doc and slides on internet about otmixmod
thanks !

9

Hi Flore,

Following our meeting here are the extended versions of the otmixmod test adapted to a multivariate input/output model and using different algorithms for the local experts:
With kriging, Q2= [0.998673,0.999992]:

import openturns as ot
import openturns.viewer as otv
import otmixmod

ot.Log.Show(ot.Log.NONE)
size = 200
kmax = 100
if size >= 1000:
    symbol = 'dot'
else:
    symbol = 'bullet'

model = ot.SymbolicFunction(
    ["x", "y", "z"], ["(1.0 + sign(x)) * cos(y) - (sign(x) - 1) * sin(2*x)", "exp(-sign(z)+x+z)"])

inDim = model.getInputDimension()
outDim = model.getOutputDimension()
# To extract the input part of the mixture
inputIndices = ot.Indices(inDim)
inputIndices.fill(0)
# To extract the output part of the mixture
outputIndices = ot.Indices(outDim)
outputIndices.fill(inDim)

dist = ot.ComposedDistribution([ot.Uniform()]*inDim)
dataX = dist.getSample(size)
dataY = model(dataX)
# For validation
sizeValid = size // 5
dataXValid = dist.getSample(sizeValid)
dataYValid = model(dataXValid)

# Stack the input and output data bases
data = ot.Sample(dataX)
data.stack(dataY)

gpairs = ot.VisualTest.DrawPairs(data)
for i in range(gpairs.getNbRows()+1):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            gpairs.setGraph(i-1, j, g)
        except:
            pass
view = otv.View(gpairs, (800, 800))
view.save("Pairs.png")
view.close()

bestExpert = 0
bestCluster = 0
bestError = ot.SpecFunc.MaxScalar
bestMixture = 0

k = 2
stop = False
covModel = 'Gaussian_pk_Lk_Ck'
if size >= 1000:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "HMAT")
else:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "LAPACK")
while not stop:
    print("Try with", k, "cluster(s)")
    try:
        # Classify data
        mixture, labels, logLike = otmixmod.MixtureFactory(k, covModel).buildAsMixture(data)
        classifier = ot.MixtureClassifier(mixture)

        # Build the clusters
        print("Build the clusters")
        clusters = otmixmod.MixtureFactory.BuildClusters(data, labels, k)

        # Build the local meta-models
        print("Build the local experts")
        metaModels = ot.Basis()
        cov = ot.MaternModel([1.0], 3.5)
        for i in range(k):
            print("Expert ", i)
            # Here we learn the kriging model output by output
            outMeta =  list()
            for n in range(outDim):
                algo = ot.KrigingAlgorithm(clusters[i].getMarginal(inputIndices),
                                           clusters[i].getMarginal(outputIndices[n]),
                                           cov, ot.Basis())
                base = ot.NLopt("LN_COBYLA")
                algo.setOptimizationAlgorithm(base)
                algo.setOptimizationBounds(ot.Interval([0.01]*inDim, [100.0]*inDim))
                algo.run()
                cov = algo.getResult().getCovarianceModel()
                outMeta.append(algo.getResult().getMetaModel())
            metaModels.add(ot.AggregatedFunction(outMeta))

        print("Build the expert mixture")
        expert = ot.Function(ot.ExpertMixture(metaModels, classifier))
        # Validation error
        dataYMeta = expert(dataXValid)
        error = (dataYValid - dataYMeta).computeRawMoment(2)
        var = dataYValid.computeCenteredMoment(2)
        for i in range(outDim):
            error[i] /= var[i]
        print("k=", k, "error=", error)
        if (error.norm() < bestError):
            bestExpert = expert
            bestCluster = k
            bestError = error.norm()
            bestMixture = mixture
    except:
        pass
    stop = (k == kmax)
    print(k, "cluster(s) done, best=", bestCluster, "error=%.4e" % bestError)
    k += 1

algo = ot.MetaModelValidation(dataXValid, dataYValid, bestExpert)
q2 = algo.computePredictivityFactor()
print("#"*50)
print("Q2=", q2)
print("#"*50)
residual = algo.getResidualDistribution()
g = algo.drawValidation()
view = otv.View(g, square_axes = True)
view.save("Validation.png")
view.close()
g = residual.drawPDF()
g.setLegendPosition("")
view = otv.View(g)
view.save("Residual.png")
view.close()
xMin = data.getMin()
xMax = data.getMax()
for i in range(1+gpairs.getNbRows()):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            g.add(bestMixture.getMarginal([j,i]).drawPDF([xMin[j], xMin[i]], [xMax[j], xMax[i]]))
            g.setLegendPosition("")
            gpairs.setGraph(i-1, j, g)
        except:
            print("failed for i=", i, "j=", j)
            pass
view = otv.View(gpairs)
view.save("PairsMixture.png")
view.close()

With basic PCE, Q2= [0.97927,0.989968]:

import openturns as ot
import openturns.viewer as otv
import otmixmod

ot.Log.Show(ot.Log.NONE)
size = 200
kmax = 100
if size >= 1000:
    symbol = 'dot'
else:
    symbol = 'bullet'

model = ot.SymbolicFunction(
    ["x", "y", "z"], ["(1.0 + sign(x)) * cos(y) - (sign(x) - 1) * sin(2*x)", "exp(-sign(z)+x+z)"])

inDim = model.getInputDimension()
outDim = model.getOutputDimension()
# To extract the input part of the mixture
inputIndices = ot.Indices(inDim)
inputIndices.fill(0)
# To extract the output part of the mixture
outputIndices = ot.Indices(outDim)
outputIndices.fill(inDim)

dist = ot.ComposedDistribution([ot.Uniform()]*inDim)
dataX = dist.getSample(size)
dataY = model(dataX)
# For validation
sizeValid = size // 5
dataXValid = dist.getSample(sizeValid)
dataYValid = model(dataXValid)

# Stack the input and output data bases
data = ot.Sample(dataX)
data.stack(dataY)

gpairs = ot.VisualTest.DrawPairs(data)
for i in range(gpairs.getNbRows()+1):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            gpairs.setGraph(i-1, j, g)
        except:
            pass
view = otv.View(gpairs, (800, 800))
view.save("Pairs.png")
view.close()

bestExpert = 0
bestCluster = 0
bestError = ot.SpecFunc.MaxScalar
bestMixture = 0

k = 2
stop = False
covModel = 'Gaussian_pk_Lk_Ck'
if size >= 1000:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "HMAT")
else:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "LAPACK")
while not stop:
    print("Try with", k, "cluster(s)")
    try:
        # Classify data
        mixture, labels, logLike = otmixmod.MixtureFactory(k, covModel).buildAsMixture(data)
        classifier = ot.MixtureClassifier(mixture)

        # Build the clusters
        print("Build the clusters")
        clusters = otmixmod.MixtureFactory.BuildClusters(data, labels, k)

        # Build the local meta-models
        print("Build the local experts")
        metaModels = ot.Basis()
        for i in range(k):
            print("Expert ", i)
            algo = ot.FunctionalChaosAlgorithm(clusters[i].getMarginal(inputIndices), clusters[i].getMarginal(outputIndices))
            algo.run()
            metaModels.add(algo.getResult().getMetaModel())
            print("distribution=", algo.getResult().getDistribution())
        print("Build the expert mixture")
        expert = ot.Function(ot.ExpertMixture(metaModels, classifier))
        # Validation error
        dataYMeta = expert(dataXValid)
        error = (dataYValid - dataYMeta).computeRawMoment(2)
        var = dataYValid.computeCenteredMoment(2)
        for i in range(outDim):
            error[i] /= var[i]
        print("k=", k, "error=", error)
        if (error.norm() < bestError):
            bestExpert = expert
            bestCluster = k
            bestError = error.norm()
            bestMixture = mixture
    except:
        pass
    stop = (k == kmax)
    print(k, "cluster(s) done, best=", bestCluster, "error=%.4e" % bestError)
    k += 1

algo = ot.MetaModelValidation(dataXValid, dataYValid, bestExpert)
q2 = algo.computePredictivityFactor()
print("#"*50)
print("Q2=", q2)
print("#"*50)
residual = algo.getResidualDistribution()
g = algo.drawValidation()
view = otv.View(g, square_axes = True)
view.save("Validation.png")
view.close()
g = residual.drawPDF()
g.setLegendPosition("")
view = otv.View(g)
view.save("Residual.png")
view.close()
xMin = data.getMin()
xMax = data.getMax()
for i in range(1+gpairs.getNbRows()):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            g.add(bestMixture.getMarginal([j,i]).drawPDF([xMin[j], xMin[i]], [xMax[j], xMax[i]]))
            g.setLegendPosition("")
            gpairs.setGraph(i-1, j, g)
        except:
            print("failed for i=", i, "j=", j)
            pass
view = otv.View(gpairs)
view.save("PairsMixture.png")
view.close()

With improved PCE, Q2= [0.997574,0.990906]:

import openturns as ot
import openturns.viewer as otv
import otmixmod

size = 200
kmax = 100
if size >= 1000:
    symbol = 'dot'
else:
    symbol = 'bullet'

model = ot.SymbolicFunction(
    ["x", "y", "z"], ["(1.0 + sign(x)) * cos(y) - (sign(x) - 1) * sin(2*x)", "exp(-sign(z)+x+z)"])

inDim = model.getInputDimension()
outDim = model.getOutputDimension()
# To extract the input part of the mixture
inputIndices = ot.Indices(inDim)
inputIndices.fill(0)
# To extract the output part of the mixture
outputIndices = ot.Indices(outDim)
outputIndices.fill(inDim)

dist = ot.ComposedDistribution([ot.Uniform()]*inDim)
dataX = dist.getSample(size)
dataY = model(dataX)
# For validation
sizeValid = size // 5
dataXValid = dist.getSample(sizeValid)
dataYValid = model(dataXValid)

# Stack the input and output data bases
data = ot.Sample(dataX)
data.stack(dataY)

gpairs = ot.VisualTest.DrawPairs(data)
for i in range(gpairs.getNbRows()+1):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            gpairs.setGraph(i-1, j, g)
        except:
            pass
view = otv.View(gpairs, (800, 800))
view.save("Pairs.png")
view.ShowAll()
view.close()

bestExpert = 0
bestCluster = 0
bestError = ot.SpecFunc.MaxScalar
bestMixture = 0

k = 2
stop = False
covModel = 'Gaussian_pk_Lk_Ck'
if size >= 1000:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "HMAT")
else:
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra", "LAPACK")
while not stop:
    print("Try with", k, "cluster(s)")
    try:
        # Classify data
        mixture, labels, logLike = otmixmod.MixtureFactory(k, covModel).buildAsMixture(data)
        classifier = ot.MixtureClassifier(mixture)

        # Build the clusters
        print("Build the clusters")
        clusters = otmixmod.MixtureFactory.BuildClusters(data, labels, k)

        # Build the local meta-models
        print("Build the local experts")
        metaModels = ot.Basis()
        for i in range(k):
            print("Expert ", i)
            # Extract the distribution of the current cluster
            distribution = mixture.getDistributionCollection()[
                i].getMarginal(inputIndices)
            # Build the local meta model using PCE
            # We use a projection strategy
            projection = ot.LeastSquaresStrategy(ot.LeastSquaresMetaModelSelectionFactory(ot.LARS(), ot.CorrectedLeaveOneOut()))
            # We use an Hermite chaos expansion
            basis = ot.OrthogonalProductPolynomialFactory(
                ot.PolynomialFamilyCollection(inDim, ot.OrthogonalUniVariatePolynomialFamily(ot.HermiteFactory())))
            # FixedStrategy
            degree = 8
            adaptive = ot.AdaptiveStrategy(
                ot.OrthogonalBasis(basis), ot.EnumerateFunction().getStrataCumulatedCardinal(degree))
            algo = ot.FunctionalChaosAlgorithm(clusters[i].getMarginal(inputIndices), clusters[
                                            i].getMarginal(outputIndices), distribution, adaptive, projection)
            algo.run()
            metaModels.add(algo.getResult().getMetaModel())

        print("Build the expert mixture")
        expert = ot.Function(ot.ExpertMixture(metaModels, classifier))
        # Validation error
        dataYMeta = expert(dataXValid)
        error = (dataYValid - dataYMeta).computeRawMoment(2)
        var = dataYValid.computeCenteredMoment(2)
        for i in range(outDim):
            error[i] /= var[i]
        print("k=", k, "error=", error)
        if (error.norm() < bestError):
            bestExpert = expert
            bestCluster = k
            bestError = error.norm()
            bestMixture = mixture
    except:
        pass
    stop = (k == kmax)
    print(k, "cluster(s) done, best=", bestCluster, "error=%.4e" % bestError)
    k += 1

algo = ot.MetaModelValidation(dataXValid, dataYValid, bestExpert)
q2 = algo.computePredictivityFactor()
print("#"*50)
print("Q2=", q2)
print("#"*50)
residual = algo.getResidualDistribution()
g = algo.drawValidation()
view = otv.View(g, square_axes = True)
view.save("Validation.png")
view.ShowAll()
view.close()
g = residual.drawPDF()
g.setLegendPosition("")
view = otv.View(g)
view.save("Residual.png")
view.ShowAll()
view.close()
xMin = data.getMin()
xMax = data.getMax()
for i in range(1+gpairs.getNbRows()):
    for j in range(i):
        try:
            g = gpairs.getGraph(i-1, j)
            d = g.getDrawable(0)
            d.setPointStyle(symbol)
            g.setDrawable(d, 0)
            g.add(bestMixture.getMarginal([j,i]).drawPDF([xMin[j], xMin[i]], [xMax[j], xMax[i]]))
            g.setLegendPosition("")
            gpairs.setGraph(i-1, j, g)
        except:
            print("failed for i=", i, "j=", j)
            pass
view = otv.View(gpairs)
view.save("PairsMixture.png")
view.ShowAll()
view.close()
10

Thanks for all Régis !!