PCE METAMODELING: CONSTRUCTION OF ORTHOGONAL BASIS FOR ARBITRARY DISTRIBUTION¶

This example demonstrates the numerical basis construction with Stieltjes procedure. The Ishigami function with Gaussian distributions restricted to the $[-\pi, \pi]$ interval is used. The Stieltjes procedure is expected to produce a basis that performs better than the Legendre polynomials, the default basis for all bounded distribution.

Package imports¶

In [1]:
from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

Initialize common plotting parameters¶

In [2]:
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()

Start a remote UQCloud session¶

In [3]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (f01a706455a94b8c89c469cf0607e892) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [4]:
uq.rng(100, 'twister');

Computational model¶

The Ishigami function is defined as:

$$Y(\mathbf{x}) = \sin(x_1) + 7 \sin^2(x_2) + 0.1 x_3^4 \sin(x_1)$$

where $x_i \in [-\pi, \pi], \; i = 1,2,3.$

This computation is carried out by the function model supplied within ishigami.py. The input parameters of this function are gathered into the vector X.

Create a MODEL from the function file:

In [5]:
ModelOpts = {
    'Type': 'Model', 
    'ModelFun':'ishigami.model'
}
In [6]:
myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The probabilistic input model consists of three independent uniform random variables:

$$X_i \sim \mathcal{U}(-\pi, \pi) \qquad i = 1,2,3 $$

Specify the marginals:

In [7]:
InputOpts = {
    "Marginals": [
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        },
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        },
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        }
    ]
}

Create an INPUT object based on the specified marginals:

In [8]:
myInput = uq.createInput(InputOpts)

PCE with Legendre polynomials¶

This section showcases several ways to calculate the polynomial chaos expansion (PCE) coefficients.

Select PCE as the metamodeling tool:

In [9]:
MetaOpts= {
    "Type": "Metamodel",
    "MetaType": "PCE",
}

Assign the Ishigami function model as the full computational model of the PCE metamodel:

In [10]:
MetaOpts["FullModel"] = myModel['Name']

The sparse-favoring least-square minimization LARS can be enabled by:

In [11]:
MetaOpts["Method"] = "LARS"

LARS allows for degree-adaptive calculation of the PCE coefficients. That is, if an array of possible degrees is given, the degree with the lowest Leave-One-Out cross-validation error (LOO error) is automatically selected. Specify the range for the degree selection:

In [12]:
MetaOpts["Degree"] = np.arange(3,16).tolist()

Manually select the Legendre orthonormal polynomials

In [13]:
MetaOpts["PolyTypes"] = ['Legendre','Legendre','Legendre']

The following options configure UQLab to generate an experimental design of size $150$ based on a latin hypercube sampling of the input model:

In [14]:
MetaOpts["ExpDesign"] = {
    "NSamples" : 150
}

Create the LARS-based PCE metamodel:

In [15]:
# reset the rng generator to obtain the same experimental design between trials
uq.rng(100, 'twister')
myPCE_Legendre = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((150,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Print a report on the calculated coefficients:

In [16]:
uq.print(myPCE_Legendre)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               15
   q-norm:                       1.00
   Size of full basis:           816
   Size of sparse basis:         50
   Full model evaluations:       150
   Leave-one-out error:          9.2344374e-12
   Modified leave-one-out error: 2.4664321e-11
   Mean value:                   3.5000
   Standard deviation:           3.7208
   Coef. of variation:           106.309%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [17]:
uq.display(myPCE_Legendre);

PCE with numerically estimated polynomials¶

We now set the polynomial families to "arbitrary", so that polynomials orthogonal to the input marginals are numerically constructed:

In [18]:
MetaOpts["PolyTypes"] = ['Arbitrary','Arbitrary','Arbitrary']

Create the OMP-based PCE metamodel:

In [19]:
# reset the rng generator to obtain the same experimental design between trials
uq.rng(100, 'twister')
myPCE_Arbitrary = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((150,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Print a report on the calculated coefficients:

In [20]:
uq.print(myPCE_Arbitrary)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               15
   q-norm:                       1.00
   Size of full basis:           816
   Size of sparse basis:         50
   Full model evaluations:       150
   Leave-one-out error:          9.2344374e-12
   Modified leave-one-out error: 2.4664321e-11
   Mean value:                   3.5000
   Standard deviation:           3.7208
   Coef. of variation:           106.309%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [21]:
uq.display(myPCE_Arbitrary);

Validation of the metamodels¶

Generation of a validation set¶

Create a validation sample of size $10^4$ from the input model

In [22]:
Xval = uq.getSample(N=1e4)

Evaluate the full model response at the validation sample points:

In [23]:
Yval = uq.evalModel(myModel,Xval)

Evaluate the corresponding responses for each of the three PCE metamodels created before:

In [24]:
YLARS = uq.evalModel(myPCE_Legendre,Xval)
YArb = uq.evalModel(myPCE_Arbitrary,Xval)
In [25]:
YPCE = [YLARS, YArb]

Comparison of the results¶

To visually assess the performance of each metamodel, produce scatter plots of the metamodel vs. the true response on the validation set:

In [26]:
methodLabels = ['Legendre','Arbitrary']

fig = plt.figure(figsize=(3,3), dpi=300, constrained_layout=True)

for idx,mLabel in enumerate(methodLabels):
    plt.subplot(1,2,idx+1)
    plt.scatter(Yval, YPCE[idx], color=uq_colors[idx], marker='.',alpha=.5, s=1)
    plt.plot([np.min(Yval), np.max(Yval)], [np.min(Yval), np.max(Yval)], 'k--',linewidth=.5)
    plt.title(methodLabels[idx],fontsize=8)
    plt.xlabel('$Y_{true}$',fontsize=6)
    plt.ylabel('$Y_{PC}$',fontsize=6)
    plt.tick_params(axis='both', labelsize=6)
    plt.xlim(np.min(Yval), np.max(Yval))
    plt.ylim(np.min(Yval), np.max(Yval))
    plt.gca().set_aspect('equal', adjustable='box')
No description has been provided for this image

Terminate the remote UQCloud session¶

In [27]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session f01a706455a94b8c89c469cf0607e892 terminated.
Out[27]:
True