Skip to content
Snippets Groups Projects
Commit c7dcaa3e authored by Nando Hegemann's avatar Nando Hegemann
Browse files

Merge branch 'add-notebooks' into 'main'

Add tutorial notebooks

See merge request !3
parents 13045566 e18d69aa
No related branches found
No related tags found
1 merge request!3Add tutorial notebooks
%% Cell type:markdown id:00c2cc0d-24ab-4a35-bc11-459feebe2455 tags:
# Tutorial 1 - Function Approximation with PyThia
%% Cell type:code id:7bd31586-38d1-4e41-8963-a0e02a070f7e tags:
``` python
import sys
import numpy as np
import pythia as pt # see doc: www.pythia-uq.com
import matplotlib.pyplot as plt
sys.path.append("../../src")
from utils import Batman # local imports
```
%% Cell type:markdown id:31a019e4-dcc1-47d9-be36-24b95e24abe9 tags:
## Define Target Function
%% Cell type:code id:ceb41891-6c8e-40ed-b55e-979e00713ea9 tags:
``` python
def target_function(y):
return Batman().eval(y)
```
%% Cell type:code id:d4a36225-6620-4919-95fa-cffe3446fef6 tags:
``` python
y = np.linspace(-7.5, 7.5, 401)
fig, ax = plt.subplot_mosaic([["batman"]], constrained_layout=True, figsize=(10, 5))
ax["batman"].plot(y, target_function(y)[:, 0], label="$f_1(y)$")
ax["batman"].plot(y, target_function(y)[:, 1], label="$f_2(y)$")
ax["batman"].set_xlabel("y")
ax["batman"].legend()
plt.show();
```
%% Cell type:markdown id:fde90e4d-be77-4b89-ae03-9ede7faa2447 tags:
## Setup Polynomial Chaos Expansion
%% Cell type:code id:b3f735c8-debd-42ba-b642-d4fb0f1ad62e tags:
``` python
param1 = pt.parameter.Parameter(index=1, name="y_1", domain=[-7.5, 7.5], distribution='uniform')
params = [param1]
```
%% Cell type:code id:67d2abd1-2ceb-442e-ae7e-417890a86e4b tags:
``` python
sdim = [75]
indices = pt.index.tensor(sdim)
index_set = pt.index.IndexSet(indices)
print(f"Multiindex information:")
print(f" number of indices: {index_set.shape[0]}")
print(f" dimension: {index_set.shape[1]}")
print(f" maximum dimension: {index_set.max}")
print(f" number of sobol indices: {len(index_set.sobol_tuples)}")
```
%% Cell type:code id:28be82de-2de1-45e6-8d83-7b61175b2e87 tags:
``` python
s = pt.sampler.ParameterSampler(params)
y_train = s.sample(1000)
w_train = s.weight(y_train)
f_train = target_function(y_train)
```
%% Cell type:markdown id:59364c66-9409-4214-8c35-4e7ed1afaa72 tags:
## Run Polynomial Chaos Approximation
%% Cell type:code id:5bea2800-edb2-4941-b238-d94825b0db23 tags:
``` python
surrogate = pt.chaos.PolynomialChaos(
params, index_set, y_train, w_train, f_train)
```
%% Cell type:markdown id:34011607-9de1-4ce0-9bd7-6b843a662703 tags:
## Compute Approximation Error
%% Cell type:code id:269991dd-bbca-4a64-a9e6-89fe2ee0f845 tags:
``` python
test_sampler = pt.sampler.ParameterSampler(params)
y_test = test_sampler.sample(1000)
f_test = target_function(y_test)
f_approx = surrogate.eval(y_test)
```
%% Cell type:code id:baa189c4-55d9-4087-8330-7de826406097 tags:
``` python
y = np.linspace(-7.1, 7.1, 200).reshape(-1, 1)
f_target = target_function(y)
f_approx = surrogate.eval(y)
f_error = np.abs(f_target-f_approx)/np.abs(f_target)
fig, ax = plt.subplot_mosaic([["target", "target"],
["approx", "error"]],
constrained_layout=True, figsize=(10, 4))
fig.suptitle("Approximation of function with PyThia")
ax["target"].scatter(y_train[:,0], f_train[:,0], s=1)
ax["target"].scatter(y_train[:,0], f_train[:,1], s=1)
ax["target"].set_title("training data points")
ax["approx"].plot(y, f_approx[:,0])
ax["approx"].plot(y, f_approx[:,1])
ax["approx"].set_title("surrogate")
ax["error"].semilogy(y, f_error[:,0])
ax["error"].semilogy(y, f_error[:,1])
ax["error"].set_title("pointwise relaive error")
plt.show()
```
%% Cell type:code id:5d1f5311-1a1e-42c7-a20e-2f8a274ce434 tags:
``` python
```
%% Cell type:markdown id:ceb7573d-54c3-417a-b533-4a87969adc59 tags:
# Tutorial 2 - Sensitivity Analysis with Sobol Indices
%% Cell type:code id:42d7cee5-4691-41c5-83d2-4ce887966690 tags:
``` python
import sys
import numpy as np
import pythia as pt # see doc: www.pythia-uq.com
sys.path.append("../../src")
from utils import Sobol # local imports
```
%% Cell type:markdown id:462aae79-a70a-4a9c-90c7-3d73b71dc287 tags:
## Define Target Function
For $y\sim\mathcal{U}([0,1]^M)$ define
$$
f(y) = \prod_{m=1}^M\frac{\vert 4y_m-2\vert + a_m}{1+a_m}
\qquad\mbox{for }a_1,\dots,a_M >= 0.
$$
%% Cell type:code id:2a2462f7-4abb-4468-89ea-a33e8268209f tags:
``` python
a = np.array([1, 2, 3, 4])
sobol = Sobol(a)
def target_function(y):
return sobol.eval(y)
```
%% Cell type:markdown id:ce631590-d241-40a9-8604-9b4035931a31 tags:
## Compute Polynomial Chaos Approximation
%% Cell type:code id:7004ecb6-18b9-4f36-96a7-b36e08a09a2d tags:
``` python
# define parameters
params = [pt.parameter.Parameter(index=j, name=f"y_{j+1}", domain=[0, 1], distribution='uniform')
for j in range(sobol.a.size)]
```
%% Cell type:code id:4911ed05-8f09-4eb9-8880-7c05e49e0f93 tags:
``` python
# set index set
indices = pt.index.simplex(len(params), 6) # limit total polynomial degree to 10
index_set = pt.index.IndexSet(indices)
print(f"multiindex information:")
print(f" number of indices: {index_set.shape[0]}")
print(f" dimension: {index_set.shape[1]}")
print(f" number of sobol indices: {len(index_set.sobol_tuples)}")
```
%% Cell type:code id:35c9b837-e3ee-436c-b7d7-c3fca55bf389 tags:
``` python
# generate training data
s = pt.sampler.ParameterSampler(params)
y_train = s.sample(10_000)
w_train = s.weight(y_train)
f_train = target_function(y_train)
```
%% Cell type:code id:212a4b5c-fff6-4d36-bb3f-c60d0304a3c2 tags:
``` python
# compute PC expansion
surrogate = pt.chaos.PolynomialChaos(params, index_set, y_train, w_train, f_train)
```
%% Cell type:code id:3c4429d6-90e6-4cf0-bb73-24e50c198f15 tags:
``` python
# test approximation
s_test = pt.sampler.ParameterSampler(params)
y_test = s_test.sample(1_000)
f_test = target_function(y_test)
f_approx = surrogate.eval(y_test)
error_L2 = np.sqrt(np.sum((f_test-f_approx)**2)/y_test.shape[0])
error_L2_rel = error_L2 / np.sqrt(np.sum((f_test)**2)/y_test.shape[0])
error_max = np.max(np.abs(f_test-f_approx))
error_max_rel = np.max(np.abs(f_test-f_approx)/np.abs(f_test))
print(f"test error L2 (abs/rel):",
f" {error_L2:4.2e} / {error_L2_rel:4.2e}")
print(f"test error max (abs/rel):",
f" {error_max:4.2e} / {error_max_rel:4.2e}")
```
%% Cell type:markdown id:bdbebf3d-d810-401f-984a-79763d8083e0 tags:
## Compare Sobol Indices
For the Sobol function we have $\mathbb{E}[f]=1$ and
$$
\operatorname{Var}[f] = \prod_{m=1}^M (1+c_m)-1
\qquad\mbox{for}\qquad
c_m = \frac{1}{3(1+a_m)^2}.
$$
The Sobol indices are thus given by
$$
\operatorname{Sob}_{m_1,\dots,m_s}
= \frac{1}{\operatorname{Var}[f]} \prod_{k=1}^s c_{m_k}.
$$
%% Cell type:code id:3c736cef-fe28-4ef6-b32a-e3da774d483d tags:
``` python
# compare Sobol indices
tuples, indices = sobol.indices()
approx = surrogate.sobol.T[0]
difference = np.abs(indices-approx)
print(f"{'tuple':<15} {'exact':<10} {'approx':<10} {'difference':<9}")
print("-"*48)
for tup, ind, app, diff in zip(tuples, indices, approx, difference):
print(f"{str(tup):<14} {ind:4.2e} {app:4.2e} {diff:4.2e}")
```
%% Cell type:code id:a99601f6-52a1-4709-94f6-0aeb2136689a tags:
``` python
```
%% Cell type:markdown id:39584f41-1647-4f93-a3b5-84a39b3e13f9 tags:
# Tutorial 3 - Markov-Chain Monte Carlo
%% Cell type:code id:914c9beb-b5de-4a26-a709-3005371211c9 tags:
``` python
import sys
import numpy as np
import pythia as pt # see doc: www.pythia-uq.com
import emcee
import matplotlib.pyplot as plt
import matplotlib.cm as cmap
sys.path.append("../../src")
from utils import Sine, Gaussian # local imports
```
%% Cell type:markdown id:43851283-9de8-42c3-90bc-1020638bdbb2 tags:
## Define Target Function
For $y\sim\mathcal{U}([0,1]^6)$ define
$$
f(y) = (f_1(y), \dots, f_{200}(y)) \in\mathbb{R}^{200}
\qquad\mbox{with}\qquad
f_j(y) = \sum_{m=1}^6 e^{-m}y_m \sin\Bigl(\frac{2\pi(j-1)}{199}\Bigr).
$$
%% Cell type:code id:31fdf19b-020e-41e4-b4ab-71c9502271c6 tags:
``` python
sine = Sine(np.exp(-np.arange(1,7)))
def target_function(y):
return sine.eval(y)
```
%% Cell type:code id:a19c32ae-3860-4194-9a13-21741a1da6db tags:
``` python
y = np.random.uniform(0, 1, (10, 6))
val = target_function(y)
fig, ax = plt.subplot_mosaic([["sine"]], constrained_layout=True, figsize=(10, 5))
for j, valj in enumerate(val):
ax["sine"].plot(sine.points, valj, color=cmap.tab20(j), label=f"sample {j+1}")
ax["sine"].legend()
plt.show();
```
%% Cell type:markdown id:c7d59cb5-3a68-4750-bf18-9220f06a224c tags:
## Set Measurement
%% Cell type:code id:9507a69e-9811-4be4-b15e-59530989018c tags:
``` python
# measurement
y_true = np.random.uniform(0.2, 0.8, (1, sine.scales.size))
f_true = target_function(y_true).reshape(1, -1)
c = 0.03
noise = np.random.normal(0, c, f_true.shape)
f_meas = f_true + noise
```
%% Cell type:markdown id:ee946804-feda-4a90-93af-e702a267c86f tags:
## Define Bayesian Posterior
%% Cell type:code id:59f218b2-2470-4672-9e69-de398956c776 tags:
``` python
domain = np.array([[0, 0, 0, 0, 0, 0, 0.01],
[1, 1, 1, 1, 1, 1, 0.5]]).T
prior = pt.sampler.ProductSampler([pt.sampler.UniformSampler(dom) for dom in domain]) # extended prior
forward_model = lambda y: target_function(y[:,:-1]) # extended forward model
sigma = lambda y: np.abs(y[:, -1]).reshape(-1, 1)*np.ones((1, f_meas.shape[1])) # error model
# define Bayesian posterior
lh = Gaussian(forward_model, sigma, xdim=domain.shape[0])
def log_posterior(y, meas):
return lh.log_likelihood(y, meas) + prior.log_pdf(y)
```
%% Cell type:markdown id:37d82fb1-9de9-4e4e-9f33-b9d10d0f2e76 tags:
## Setup MCMC
%% Cell type:code id:9e867a60-7337-4f24-ab14-e28f60285376 tags:
``` python
n_walkers = 15
chain_length = 2000
seed = prior.sample(n_walkers)
```
%% Cell type:code id:fa098095-592e-4b76-ba86-8515fb9b0c76 tags:
``` python
# define runner
runner = emcee.EnsembleSampler(n_walkers, domain.shape[0], log_posterior, args=[f_meas])
```
%% Cell type:markdown id:76308649-ecad-4a4c-ade5-8074f7673e32 tags:
## Run MCMC
%% Cell type:code id:06aedd14-aadf-4f14-a6b3-3a32c0286286 tags:
``` python
# discard burn-in
state = runner.run_mcmc(seed, 1000, progress=True)
runner.reset()
```
%% Cell type:code id:2ab5096c-8e4a-4dc3-8edd-02a543c495ce tags:
``` python
# run MCMC
runner.run_mcmc(seed, chain_length, progress=True);
```
%% Cell type:code id:5b2684a6-1a92-4820-a5da-d9259a074d49 tags:
``` python
# simple evaluation
samples = runner.get_chain(flat=True)
mean = np.mean(samples, axis=0)
std = np.std(samples, axis=0)
```
%% Cell type:markdown id:14dfbf69-c251-43dd-bfa2-5588b0bcc3cc tags:
## Data Analysis
%% Cell type:code id:b84800c1-350a-4581-975c-58e238c83b91 tags:
``` python
# plot posterior distribution
layout = [[0, 1], [2, 3], [4, 5], [6, 6]]
fig, axes = plt.subplot_mosaic(layout, figsize=(12, 7.5))
for j, (dom, ax) in enumerate(zip(domain, axes.values())):
bins = 100
ax.hist(samples[:, j], bins=bins, density=True, range=dom, alpha=0.5)
ax.axvline(x=mean[j], color=cmap.tab10(0), label="mean")
ax.axvline(x=mean[j]-std[j], ls="--", color=cmap.tab10(0), label="std")
ax.axvline(x=mean[j]+std[j], ls="--", color=cmap.tab10(0))
if j == 6:
ax.axvline(x=c, color="k", label="true value")
ax.set_ylabel(f"$c$")
else:
ax.axvline(x=y_true[0, j], color="k", label="true value")
ax.set_ylabel(f"$y_{j+1}$")
ax.set_yticks([])
ax.legend()
plt.show()
```
%% Cell type:code id:9f3e4661-663e-41ab-beb4-5734a11a4d94 tags:
``` python
```
%% Cell type:markdown id:b304f6f3-ad5b-4123-9560-316b6d5cb434 tags:
# Tutorial 4 - Sampling
%% Cell type:code id:78eb2cc1-33e6-4fe3-bb47-70e6aaca324a tags:
``` python
import numpy as np
import pythia as pt # check doc: www.pythia-uq.com
import matplotlib.pyplot as plt
import matplotlib.cm as cmap
```
%% Cell type:markdown id:9b31d394-d268-4c7e-af46-bff7c38d71e2 tags:
## Univariate Sampling
%% Cell type:code id:69a027e1-e887-47ea-882b-04d2781b4a93 tags:
``` python
# uniform
s = pt.sampler.UniformSampler([-1,2])
samples = s.sample(100_000)
y = np.linspace(-1.5, 2.5, 200)
fig, ax = plt.subplot_mosaic([["fig"]], constrained_layout=True)
ax["fig"].hist(samples.flatten(), bins=50, density=True, label="samples")
ax["fig"].plot(y, s.pdf(y), label="pdf")
ax["fig"].legend()
plt.show()
```
%% Cell type:code id:bda40e48-deb9-4095-a25c-f5fd3b7dd5b5 tags:
``` python
# normal
s = pt.sampler.NormalSampler(mean=1, var=1)
samples = s.sample(100_000)
y = np.linspace(-3.5, 5.5, 200)
fig, ax = plt.subplot_mosaic([["fig"]], constrained_layout=True)
ax["fig"].hist(samples.flatten(), bins=50, density=True, label="samples")
ax["fig"].plot(y, s.pdf(y), label="pdf")
ax["fig"].legend()
plt.show()
```
%% Cell type:code id:c0dea9d0-1c87-43c7-be74-cebba7efe6c1 tags:
``` python
# Gamma
s1 = pt.sampler.GammaSampler(domain=[0, np.inf], alpha=1, beta=1)
samples1 = s1.sample(100_000)
s2 = pt.sampler.GammaSampler(domain=[0, np.inf], alpha=5, beta=5)
samples2 = s2.sample(100_000)
y = np.linspace(-0.5, 5.5, 200)
fig, ax = plt.subplot_mosaic([["a", "b"]], constrained_layout=True)
ax["a"].hist(samples1.flatten(), bins=50, density=True, label="samples")
ax["a"].plot(y, s1.pdf(y), label="pdf")
ax["a"].legend()
ax["b"].hist(samples2.flatten(), bins=50, density=True, label="samples")
ax["b"].plot(y, s2.pdf(y), label="pdf")
ax["b"].legend()
plt.show()
```
%% Cell type:code id:8e9ab462-81d0-4f9c-9fa9-21884ba1a9b1 tags:
``` python
# Beta
s1 = pt.sampler.BetaSampler(domain=[0, 3], alpha=0.8, beta=0.8)
samples1 = s1.sample(100_000)
s2 = pt.sampler.BetaSampler(domain=[0, 3], alpha=10, beta=2)
samples2 = s2.sample(100_000)
y = np.linspace(-0.5, 3.5, 200)
fig, ax = plt.subplot_mosaic([["a", "b"]], constrained_layout=True)
ax["a"].hist(samples1.flatten(), bins=50, density=True, label="samples")
ax["a"].plot(y, s1.pdf(y), label="pdf")
ax["a"].legend()
ax["b"].hist(samples2.flatten(), bins=50, density=True, label="samples")
ax["b"].plot(y, s2.pdf(y), label="pdf")
ax["b"].legend()
plt.show()
```
%% Cell type:markdown id:69e22499-9770-47e1-becb-326e29746947 tags:
## Multivariate Sampling
%% Cell type:code id:9c2a8fb8-239f-40a0-8fc4-d9776eb1c4a8 tags:
``` python
# product sampler
s1 = pt.sampler.UniformSampler(domain=[0, 1])
s2 = pt.sampler.NormalSampler(mean=0, var=1)
s = pt.sampler.ProductSampler([s1, s2])
samples = s.sample(100_000)
x, y = np.linspace(0, 1, 50), np.linspace(-4, 4, 50)
xy = pt.misc.cartProd([x, y])
pdf = s.pdf(xy).reshape(x.size, y.size).T
fig, ax = plt.subplot_mosaic([["x1", "pdf"], ["samples", "x2"]], constrained_layout=True)
ax["pdf"].contourf(pdf, 30, cmap="Greys", extent=[0, 1, -4, 4])
ax["samples"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color="grey")
ax["samples"].axhline(y=0, c=cmap.tab10(0))
ax["samples"].axvline(x=0.5, c=cmap.tab10(1))
ax["x1"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))
ax["x2"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))
plt.show();
```
%% Cell type:code id:7903b7f1-36a4-4899-a3b2-2ffbc8abce8b tags:
``` python
# parameter sampler
param1 = pt.parameter.Parameter(index=1, name="y1", domain=[0, 1], distribution="uniform")
param2 = pt.parameter.Parameter(index=2, name="y2", domain=[-np.inf, np.inf], distribution="normal",
mean=0, var=1)
params = [param1, param2]
s = pt.sampler.ParameterSampler(params)
samples = s.sample(100_000)
x, y = np.linspace(0, 1, 50), np.linspace(-4, 4, 50)
xy = pt.misc.cartProd([x, y])
pdf = s.pdf(xy).reshape(x.size, y.size).T
fig, ax = plt.subplot_mosaic([["x1", "pdf"], ["samples", "x2"]], constrained_layout=True)
ax["pdf"].contourf(pdf, 30, cmap="Greys", extent=[0, 1, -4, 4])
ax["samples"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color="grey")
ax["samples"].axhline(y=0, c=cmap.tab10(0))
ax["samples"].axvline(x=0.5, c=cmap.tab10(1))
ax["x1"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))
ax["x2"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))
plt.show();
```
%% Cell type:code id:ba9900d0-b2e3-4bf0-9bbb-16c09df1eafc tags:
``` python
# WLS tensor sampler
param1 = pt.parameter.Parameter(index=1, name="y1", domain=[0, 1], distribution="uniform")
param2 = pt.parameter.Parameter(index=2, name="y2", domain=[0, 1], distribution="beta",
alpha=5, beta=5)
params = [param1, param2]
s = pt.sampler.WLSTensorSampler(params, [3, 3], tsa=True)
samples = s.sample(100_000)
x, y = np.linspace(0, 1, 50), np.linspace(0, 1, 50)
xy = pt.misc.cartProd([x, y])
pdf = s.pdf(xy).reshape(x.size, y.size).T
fig, ax = plt.subplot_mosaic([["x1", "pdf"], ["samples", "x2"]], constrained_layout=True)
ax["pdf"].contourf(pdf, 30, cmap="Greys", extent=[0, 1, 0, 1])
ax["samples"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color="grey")
ax["samples"].axhline(y=0.5, c=cmap.tab10(0))
ax["samples"].axvline(x=0.5, c=cmap.tab10(1))
ax["x1"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))
ax["x2"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))
plt.show();
```
%% Cell type:code id:1dc2824d-e495-453b-9eba-f48bfdc03945 tags:
``` python
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment