Get the output marginal of a PCE

pce_marginal.py

"""
Get the (output) marginal of a PCE.

There is no FunctionalChaosResult.getMarginal()
https://github.com/openturns/openturns/issues/2985

Limitation
----------

This script involves more PCE indices than necessary i.e. the number 
of coefficients matches one for all output marginals. 
This is because it does not remove the indices for which the 
marginal coefficient is zero.

Example
-------

We consider the flood model:
- inputs : [Q,Ks,Zv,Zm,B,L,Zb,Hd]
- output : [H,S,C]

We create a PCE for all outputs. 
Then we want to get the 3 marginal PCEs for H, S and C.

The next script prints:

Marginal PCE #0 / 3
FunctionalChaosResult
- input dimension=8
- output dimension=1
- distribution dimension=8
- transformation=8 -> 8
- inverse transformation=8 -> 8
- orthogonal basis dimension=8
- indices size=28
- relative errors=[0.000208468]
- residuals=[0.00699816]
- is least squares=true
- is model selection=false

| Index | Multi-index   | Coefficient   |
|-------|---------------|---------------|
|     0 | [0,0,0,0,0,0,0,0]| 2.53894       |
[...]
|    27 | [0,0,0,1,7,0,0,0]| 0             |

Marginal PCE #1 / 3
FunctionalChaosResult
- input dimension=8
[...]
- is model selection=false

| Index | Multi-index   | Coefficient   |
|-------|---------------|---------------|
|     0 | [0,0,0,0,0,0,0,0]| -5.89789      |
[...]
|    27 | [0,0,0,1,7,0,0,0]| 0             |

Marginal PCE #2 / 3
FunctionalChaosResult
- input dimension=8
- output dimension=1
[...]
- is model selection=false

| Index | Multi-index   | Coefficient   |
|-------|---------------|---------------|
|     0 | [0,0,0,0,0,0,0,0]| 1.04658       |
[...]
|    27 | [0,0,0,1,7,0,0,0]| 0.0116759     |
"""

# %%
from openturns.usecases import flood_model
import openturns as ot

# %%
fm = flood_model.FloodModel()
sampleSize = 100
inputTrain = fm.distribution.getSample(sampleSize)
print(f"Output dimension = {fm.model.getOutputDimension()}")
outputTrain = fm.model(inputTrain)
marginalList = [fm.distribution.getMarginal(i) for i in range(fm.distribution.getDimension())]
multivariateBasis = ot.OrthogonalProductPolynomialFactory(marginalList)
selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory()
projectionStrategy = ot.LeastSquaresStrategy(selectionAlgorithm)
totalDegree = 8
enumerateFunction = multivariateBasis.getEnumerateFunction()
basisSize = enumerateFunction.getBasisSizeFromTotalDegree(totalDegree)
adaptiveStrategy = ot.FixedStrategy(multivariateBasis, basisSize)
chaosAlgo = ot.FunctionalChaosAlgorithm(
    inputTrain, outputTrain, fm.distribution, adaptiveStrategy, projectionStrategy
)
chaosAlgo.run()
chaosResult = chaosAlgo.getResult()
chaosResult


# %%
# Extract PCE marginal
def getMarginalPCE(chaosResult, index_output=0):
    """Extract the marginal Polynomial Chaos Expansion (PCE) for a specific output.

    Parameters
    ----------
    chaosResult : ot.FunctionalChaosResult
        Result from a Functional Chaos analysis.
    index_output : int
        Index of the output component to extract.

    Returns
    -------
    marginal_pce : ot.FunctionalChaosResult
        Marginal PCE for the specified output component.
    """
    input_sample = chaosResult.getInputSample()  # 1
    output_sample = chaosResult.getOutputSample()  # 2
    distribution = chaosResult.getDistribution()  # 3
    transformation = chaosResult.getTransformation()  # 4
    inverse_transformation = chaosResult.getInverseTransformation()  # 5
    basis = chaosResult.getOrthogonalBasis()  # 6
    indices = chaosResult.getIndices()  # 7
    coefficients = chaosResult.getCoefficients()  # 8
    functions = chaosResult.getReducedBasis()  # 9
    residuals_point = chaosResult.getResiduals()  # 10
    relative_errors_point = chaosResult.getRelativeErrors()  # 11

    marginal_pce = ot.FunctionalChaosResult(
        input_sample, # 1
        output_sample.getMarginal(index_output), # 2
        distribution, # 3
        transformation, # 4
        inverse_transformation, # 5
        basis, # 6
        indices, # 7 TODO: simplify this to get the smallest possible indices
        coefficients.getMarginal(index_output), # 8
        functions, # 9
        [residuals_point[index_output]], # 10
        [relative_errors_point[index_output]] # 11
    )
    return marginal_pce

# %%
marginal_pce_results = []
for i in range(fm.model.getOutputDimension()):
    marginal_pce = getMarginalPCE(chaosResult, i)
    marginal_pce_results.append(marginal_pce)

# %%
for i in range(fm.model.getOutputDimension()):
    print(f"Marginal PCE #{i} / {fm.model.getOutputDimension()}:")
    print(marginal_pce_results[i])

# %%