PCE METAMODELING: EXPERIMENTAL DESIGN OPTIONS¶

This example shows different methods to create an experimental design and calculate PCE coefficients with the sparse-favoring least-square minimization strategy. The full computational model of choice is the Ishigami function, a standard benchmark in PCE applications.

Package imports¶

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

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 (be9781fafac942c5a1371df02938770e) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [4]:
uq.rng(120,'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 $N \times M$ matrix X, where $N$ and $M$ are the numbers of realizations and input variables, respectively.

Create a MODEL from the function file:

In [5]:
ModelOpts = {
    'Type': 'Model', 
    'ModelFun':'ishigami.model'
}
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 [6]:
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 [7]:
myInput = uq.createInput(InputOpts)

Polynomial chaos expansion (PCE) metamodels¶

In this section, a sparse polynomial chaos expansion (PCE) is created to surrogate the Ishigami function. The method of choice is the sparse-favouring least-square minimization strategy LARS.

Select PCE as the metamodeling tool and LARS as the calculation strategy:

In [8]:
MetaOpts = {
    "Type": "Metamodel",
    "MetaType": "PCE",
    "Method": "LARS"
}

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

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

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 (LOO) cross-validation error is automatically selected. Specify the range for the degree selection:

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

Comparison of the performance of the metamodels created below is done by comparing the response of the full model and each of the metamodels on a validation set. Generate a validation set of size $10^4$:

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

Calculate the response of the exact Ishigami model at the validation set points for reference:

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

COMPARISON OF DIFFERENT EXPERIMENTAL DESIGN SIZES¶

Least-square methods rely on the evaluation of the model response on an experimental design. The following options configure UQLab to generate an experimental design of size $20$ based on a latin hypercube sampling of he input model (also available: 'MC', 'Sobol', 'Halton'):

In [13]:
MetaOpts["ExpDesign"] = {
    "NSamples":20,
    "Sampling":"LHS"
}

Create the PCE metamodel based on the design of size $20$:

In [14]:
myPCE_LHS_20 = 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 ((20,))...

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

Create additional PCE metamodels with an increasing number of points in the experimental design: $40$, $80$, and $120$ sample points; with all the other options are left unchanged:

In [15]:
MetaOpts["ExpDesign"]["NSamples"] = 40
myPCE_LHS_40 = uq.createModel(MetaOpts)

MetaOpts["ExpDesign"]["NSamples"] = 80
myPCE_LHS_80 = uq.createModel(MetaOpts)

MetaOpts["ExpDesign"]["NSamples"] = 120
myPCE_LHS_120 = 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 ((40,))...

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

 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 ((80,))...

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

 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 ((120,))...

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

Calculate the responses of the created metamodels:

In [16]:
Y_LHS_20 = uq.evalModel(myPCE_LHS_20,Xval)
Y_LHS_40 = uq.evalModel(myPCE_LHS_40,Xval)
Y_LHS_80 = uq.evalModel(myPCE_LHS_80,Xval)
Y_LHS_120 = uq.evalModel(myPCE_LHS_120,Xval)
Y_LHS = [Y_LHS_20, Y_LHS_40, Y_LHS_80, Y_LHS_120]

Create plots for visual comparison between all the metamodels¶

created with different sample sizes:

In [17]:
methodLabels = ["N=20", "N=40", "N=80", "N=120"]

for idx,mLabel in enumerate(methodLabels):
    plt.subplot(2,2,idx+1)
    plt.scatter(Yval, Y_LHS[idx], marker='+',alpha=.5)
    plt.plot([np.min(Yval), np.max(Yval)], [np.min(Yval), np.max(Yval)], 'k')
    plt.title(methodLabels[idx])
    plt.grid('True')
    if idx > 1  : plt.xlabel('$Y_{true}$')
    if idx % 2 == 0 : plt.ylabel('$Y_{PC}$')
    plt.tick_params(axis='both')
    plt.xlim(np.min(Yval), np.max(Yval))
    plt.ylim(np.min(Yval), np.max(Yval))
    plt.gca().set_aspect('equal', adjustable='box')
plt.tight_layout()
plt.show()
No description has been provided for this image

COMPARISON OF DIFFERENT SAMPLING STRATEGIES¶

Create additional metamodels with experimental designs of fixed size ($N = 80$) based on different sampling strategies (Monte Carlo sampling,Sobol', and Halton' sequences):

In [18]:
MetaOpts["ExpDesign"]["NSamples"] = 80

MetaOpts["ExpDesign"]["Sampling"] = 'MC'
myPCE_MC_80 = uq.createModel(MetaOpts)

MetaOpts["ExpDesign"]["Sampling"] = 'Sobol'
myPCE_Sobol_80 = uq.createModel(MetaOpts)

MetaOpts["ExpDesign"]["Sampling"] = 'Halton'
myPCE_Halton_80 = 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 ((80,))...

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

 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 ((80,))...

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

 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 ((80,))...

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

Calculate the responses of the created metamodels:

In [19]:
Y_MC_80 = uq.evalModel(myPCE_MC_80,Xval)
Y_Sobol_80 = uq.evalModel(myPCE_Sobol_80,Xval)
Y_Halton_80 = uq.evalModel(myPCE_Halton_80,Xval)
Y_80 = [Y_LHS_80, Y_MC_80, Y_Sobol_80, Y_Halton_80]

Create plots for visual comparison between all the metamodels¶

created with different sampling strategies:

In [20]:
methodLabels = ["LHS", "MC", "Sobol'", "Halton"]

for idx,mLabel in enumerate(methodLabels):
    plt.subplot(2,2,idx+1)
    plt.scatter(Yval, Y_80[idx], marker='+',alpha=.5)
    plt.plot([np.min(Yval), np.max(Yval)], [np.min(Yval), np.max(Yval)], 'k')
    plt.title(methodLabels[idx])
    plt.grid('True')
    if idx > 1  : plt.xlabel('$Y_{true}$')
    if idx % 2 == 0 : plt.ylabel('$Y_{PC}$')
    plt.tick_params(axis='both')
    plt.xlim(np.min(Yval), np.max(Yval))
    plt.ylim(np.min(Yval), np.max(Yval))
    plt.gca().set_aspect('equal', adjustable='box')
plt.tight_layout()
plt.show()
No description has been provided for this image

Terminate the remote UQCloud session¶

In [21]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session be9781fafac942c5a1371df02938770e terminated.
Out[21]:
True