Tutorial - Shape Reconstruction via Scatterometry

This tutorial shows how PyThia can be used to approximate the scatterometry measurement process and conduct a global sensitivity analysis of the model.

In this experiment a silicon line grating is irradiated by an electromagnetic wave under different angles of incidence \(\theta\) and azimuthal orientations \(\phi\). The aim of this experiment is to build a surrogate for the map from the six shape parameters \(x = (\mathit{h}, \mathit{cd}, \mathit{swa}, \mathit{t}, \mathit{r}_{\mathrm{top}}, \mathit{r}_{\mathrm{bot}})\) to the detected scattered far field of the electromagnetic wave. The experimental setup is visualized in the image below.



If you want to run the source code for this tutorial yourself, you need to download the tutorial data. This tutorial assumes that the data are located in the same directory as the source code.


We start by loading the numpy and pythia packages.

import numpy as np
import pythia as pt

Generating training data from the forward problem requires solving Maxwell’s equations with finite elements, which is very intricate and time consuming. We assume the training data we rely on in this tutorial have been previously generated in an “offline” step, on a different machine using the JCMsuite software. Hence we can simply load the data.

x_data = np.load("x_train.npy")  # (10_000, 6)
y_data = np.load("y_train.npy")  # (90, 10_000, 2)
w_data = np.load("w_train.npy")  # (10_000, )

For the simulations we assumed Beta distributed parameters in reasonable parameter domains, which is also the density we drew the samples from. In particular, this implies that the weights are simply \(w_i = 1 / N\), where \(N\) is the number of samples (i.e., \(10.000\) in this case). The scattered intensity signal data consist of 4 curves (2 different azimuthal angles \(\phi\) with s- and p- polarized wave signals each) where each curve is evaluated in 45 points (angles of incidence \(\theta\)). The continuous curves look similar to the image below.


As a first step to compute a polynomial chaos expansion, we need to define the parameters using PyThias Parameter class. We can use the parameter.npy file to access the parameter information.

params = []
for dct in np.load("parameter.npy", allow_pickle=True):

Below you can find a table and an image listing the relevant parameter information.

Parameter Information








\(\mathrm{Beta}_{\alpha,\beta}([43, 53])\)





\(\mathrm{Beta}_{\alpha,\beta}([22, 28])\)





\(\mathrm{Beta}_{\alpha,\beta}([84, 90])\)





\(\mathrm{Beta}_{\alpha,\beta}([4, 6])\)





\(\mathrm{Beta}_{\alpha,\beta}([8, 13])\)





\(\mathrm{Beta}_{\alpha,\beta}([3, 7])\)




Next, we split the available data into a training and a test set, using roughly \(8.000\) data points for the training.

# split training/test dataset
split = int(0.8 * x_data.shape[0])

# training data
x_train = x_data[:split]
y_train = y_data[:split]
w_train = w_data[:split]

# test data
x_test = x_data[split:]
y_test = y_data[split:]

Now that we have our parameters defined and the training data in place, we only need to define which expansion terms we want to include into our expansion. In this case, we choose a simplex index set to bound the total polynomial degree of the expansion by four. This leaves us with \(210\) polynomial chaos terms in our expansion.

indices = pt.index.simplex_set(dimension=len(params), maximum=4)
index_set = pt.index.IndexSet(indices)

Finally, we have all the ingredients to compute our polynomial chaos surrogate, which requires only one line of code.

surrogate = pt.chaos.PolynomialChaos(params, index_set, x_train, w_train, y_train)

Evaluating the surrogate on our test set with

y_approx = surrogate.eval(x_test)

allows us compute an empirical approximation of the \(L^2\)-error (MSE), which is approximately at \(5*10^{-3}\). This is ok, since the measurement error of the observation we will investigate for the inverse problem will be at least one order of magnitude larger.

Sensitivity analysis

There are essentially two different types of sensitivity analyses, local and global. Local sensitivity analysis typicalls considers a change in the signal by a small variation in the input around a specific point. Given the polynomial chaos surrogate, computation of the local sensitivities is of course very easy as derivatives of the polynomial are given analytically. To compute the local sensitivities of the second parameter (\(\mathit{cd}\)) in all points of our test set, i.e. \(\partial/\partial x_{\mathit{cd}}\ f(x_{\mathrm{test}})\), we can use the optional partial argument of the eval() function

local_sensitivity = surrogate.eval(x_test, partial=[0, 1, 0, 0, 0, 0])

However, computing local sensitivities yields only a very limited image of the signal dependency on change in the individual parameter or combinations thereof. When it comes to estimating whether the reconstruction of the hidden parameters will be possible it is often more reliable to consider the overall dependence of the signal on the parameters. There are various types of global sensitivity analysis methods, but essentially all of them rely on integrating the target function over the high-dimensional parameter space, which is typically not feasible. The polynomial chaos expansion, however, is closely connected to the so called Sobol indices, which allows us to compute these indices by squaring and summing different subsets of our expansion coefficients. The Sobol coefficients are automatically computed with the polynomial chaos approximation and can be accessed via the surrogate.sobol attribute, which is ordered according to the list of Sobol tuples in index_set.sobol_tuples. As the signal data consist of multiple points on different curves, we get the Sobol indices for each point of our signal. As this is hard to visualize however, we can have a look on the maximal value each Sobol index takes over the signal points. The 10 most significant Sobol indices are listed in the table below.

Most relevant Sobol indices


Sobol tuple

Parameter combinations

Sobol index value (max)


















\((2, 5)\)

\((\mathit{cd}, \mathit{r}_{\mathrm{top}})\)







\((3, 5)\)

\((\mathit{swa}, \mathit{r}_{\mathrm{top}})\)



\((1, 2)\)

\((\mathit{h}, \mathit{cd})\)







\((2, 3, 5)\)

\((\mathit{cd}, \mathit{swa}, \mathit{r}_{\mathrm{top}})\)


We see that all single parameter Sobol indices are among the top ten, implying that changes in single parameters are among the most influential. Moreover, looking at the actual (maximal) values, we can assume that the first 4 parameters (\(\mathit{h}\), \(\mathit{cd}\), \(\mathit{swa}\), \(\mathit{t}\)) can be reconstructed very well, whereas the sensitivities of the corner roundings \(\mathit{r}_{\mathrm{top}}\) and \(\mathit{r}_{\mathrm{bot}}\) imply that we will get large uncertainties in the reconstruction.

This concludes the tutorial.

Complete Script

The complete source code listed below is located in the tutorials/ subdirectory of the repository root.

  2File: tutorials/tutorial_scat.py
  3Author: Nando Hegemann
  4Gitlab: https://gitlab.com/Nando-Hegemann
  5Description: Tutorial - Shape Reconstruction via Scatterometry.
  6SPDX-License-Identifier: LGPL-3.0-or-later OR Hippocratic-3.0-ECO-MEDIA-MIL
  8import numpy as np
  9import pythia as pt
 11angles = np.load("angles.npy")  # (90, 2)
 13# reshape
 14angles = np.concatenate([angles, angles], axis=0)  # (180, 2)
 15angle_idxs = [range(45), range(45, 90), range(90, 135), range(135, 180)]
 16angle_labels = [
 17    "$\\phi=0^\\circ$, s pol",
 18    "$\\phi=0^\\circ$, p pol",
 19    "$\\phi=90^\\circ$, s pol",
 20    "$\\phi=90^\\circ$, p pol",
 22print("angle info:")
 23print(f"    azimuth phi: {np.unique(angles[:,0])}")
 24print(f"    incidence theta: {np.unique(angles[:,1])}")
 26# load data
 27x_data = np.load("x_train.npy")  # (10_000, 6)
 28y_data = np.load("y_train.npy")  # (90, 10_000, 2)
 29w_data = np.load("w_train.npy")  # (10_000, )
 31# reshape y_data
 32y_data = np.concatenate(
 33    [y_data[:45, :, 0], y_data[:45, :, 1], y_data[45:, :, 0], y_data[45:, :, 1]], axis=0
 34).T  # (10_000, 180)
 35print("data info:")
 36print(f"    x_data.shape: {x_data.shape}")
 37print(f"    w_data.shape: {w_data.shape}")
 38print(f"    y_data.shape: {y_data.shape}")
 40# load parameters
 41params = []
 42for dct in np.load("parameter.npy", allow_pickle=True):
 43    params.append(
 44        pt.parameter.Parameter(
 45            name=dct["_name"],
 46            domain=dct["_dom"],
 47            distribution=dct["_dist"],
 48            alpha=dct["_alpha"],
 49            beta=dct["_beta"],
 50        )
 51    )
 53print("parameter info:")
 54print("    parameter  dist    domain    alpha   beta")
 55print("    " + "-" * 42)
 56for param in params:
 57    print(
 58        f"    {param.name:<8}  ",
 59        f"{param.distribution:<5}  ",
 60        f"{str(np.array(param.domain, dtype=int)):<7}  ",
 61        f"{param.alpha:<4.2f}  ",
 62        f"{param.beta:<4.2f}",
 63    )
 65# split training/test dataset
 66split = int(0.8 * x_data.shape[0])
 68# training data
 69x_train = x_data[:split]
 70y_train = y_data[:split]
 71w_train = w_data[:split]
 73# test data
 74x_test = x_data[split:]
 75y_test = y_data[split:]
 77# define multiindex set
 78indices = pt.index.simplex_set(dimension=len(params), maximum=4)
 79index_set = pt.index.IndexSet(indices)
 81print("Multiindex information:")
 82print(f"    number of indices: {index_set.shape[0]}")
 83print(f"    dimension: {index_set.shape[1]}")
 84print(f"    maximum dimension: {index_set.max}")
 85print(f"    number of sobol indices: {len(index_set.sobol_tuples)}")
 87print("compute PC surrogate")
 88# run pc approximation
 89surrogate = pt.chaos.PolynomialChaos(params, index_set, x_train, w_train, y_train)
 90# test approximation error
 91y_approx = surrogate.eval(x_test)
 93# L2-L2 error
 94err_abs = np.sum(np.linalg.norm(y_test - y_approx, axis=1)) / y_test.shape[0]
 95err_rel = np.sum(np.linalg.norm((y_test - y_approx) / y_test, axis=1)) / y_test.shape[0]
 96print(f"    L2-L2 error: {err_abs:4.2e} (abs), {err_rel:4.2e} (rel)")
 98# C0-L2 error
 99err_abs = np.sum(np.max(np.abs(y_test - y_approx), axis=1)) / y_test.shape[0]
100err_rel = (
101    np.sum(np.max(np.abs(y_test - y_approx) / np.abs(y_test), axis=1)) / y_test.shape[0]
103print(f"    L2-C0 error: {err_abs:4.2e} (abs), {err_rel:4.2e} (rel)")
105# global sensitivity analysis
106print("maximal Sobol indices:")
107max_vals = np.max(surrogate.sobol, axis=1)
108l2_vals = np.linalg.norm(surrogate.sobol, axis=1)
109sobol_max = [
110    list(reversed([x for _, x in sorted(zip(max_vals, index_set.sobol_tuples))]))
112sobol_L2 = [
113    list(reversed([x for _, x in sorted(zip(l2_vals, index_set.sobol_tuples))]))
115print(f"    {'#':>2}  {'max':<10}  {'L2':<10}")
116print("    " + "-" * 25)
117for j in range(10):
118    print(f"    {j+1:>2}  {str(sobol_max[0][j]):<10}  {str(sobol_L2[0][j]):<10}")