From da2f743a0911336a4e30f8337e30e45c790703ec Mon Sep 17 00:00:00 2001 From: Hu Zhao <zhao@mbd.rwth-aachen.de> Date: Thu, 8 Dec 2022 17:00:21 +0100 Subject: [PATCH] docs: implement example for ScalarGaSP --- docs/examples/emulator/README.rst | 2 - docs/examples/emulator/robustgasp/README.rst | 2 + .../emulator/robustgasp/plot_scalargasp.py | 107 ++++++++++++++++++ docs/source/examples.rst | 3 +- 4 files changed, 111 insertions(+), 3 deletions(-) delete mode 100644 docs/examples/emulator/README.rst create mode 100644 docs/examples/emulator/robustgasp/README.rst create mode 100644 docs/examples/emulator/robustgasp/plot_scalargasp.py diff --git a/docs/examples/emulator/README.rst b/docs/examples/emulator/README.rst deleted file mode 100644 index 2665689..0000000 --- a/docs/examples/emulator/README.rst +++ /dev/null @@ -1,2 +0,0 @@ -emulator Examples -^^^^^^^^^^^^^^^^^ \ No newline at end of file diff --git a/docs/examples/emulator/robustgasp/README.rst b/docs/examples/emulator/robustgasp/README.rst new file mode 100644 index 0000000..20c0480 --- /dev/null +++ b/docs/examples/emulator/robustgasp/README.rst @@ -0,0 +1,2 @@ +Emulator - RobustGaSP +--------------------- \ No newline at end of file diff --git a/docs/examples/emulator/robustgasp/plot_scalargasp.py b/docs/examples/emulator/robustgasp/plot_scalargasp.py new file mode 100644 index 0000000..e4143bc --- /dev/null +++ b/docs/examples/emulator/robustgasp/plot_scalargasp.py @@ -0,0 +1,107 @@ +""" + +ScalarGaSP: GP emulation for a single-output function +===================================================== +""" + +# %% md +# +# This example shows how apply Gaussian process emulation to a single-output +# function using class :class:`.ScalarGaSP`. +# +# The task is to build a GP emulator for the function :math:`y = x * sin(x)` +# based on a few number of training data. + +# %% md +# +# First, import the class :class:`.ScalarGaSP` by + +from psimpy.emulator import ScalarGaSP + +# %% md +# +# Then, create an instance of :class:`.ScalarGaSP`. The parameter ``ndim`` +# (dimension of function input ``x``) must be specified. Optional parameters, such +# as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave +# all the optional parameters to their default values. + +emulator = ScalarGaSP(ndim=1) + +# %% md +# +# Given training input points ``design`` and corresponding output values ``response``, +# the emulator can be trained. Below we train an emulator using :math:`8` selected +# points. + +import numpy as np + +def f(x): + #return x + 3*np.sin(x/2) + return x*np.sin(x) + +x = np.arange(2,10,1) +y = f(x) + +emulator.train(design=x, response=y) + +# %% md +# +# We can validate the performance of the trained emulator using the leave-one-out +# cross validation method :py:meth:`loo_validate()`. + +validation = emulator.loo_validate() + +# %% md +# +# Let's plot emulator predictions vs actual outputs. The error bar indicates the +# 95% credible interval. + +import matplotlib.pyplot as plt + +fig , ax = plt.subplots(figsize=(4,4)) + +ax.set_xlabel('Actual y') +ax.set_ylabel('Emulator-predicted y') +ax.set_xlim(np.min(y)-1,np.max(y)+1) +ax.set_ylim(np.min(y)-1,np.max(y)+1) + +_ = ax.plot([np.min(y)-1,np.max(y)+1], [np.min(y)-1,np.max(y)+1]) +_ = ax.errorbar(y, validation[:,0], validation[:,1], fmt='.', linestyle='') + + +# %% md +# +# With the trained emulator at our deposit, we can use the +# :py:meth:`.ScalarGaSP.predict()` method to make predictions at +# any arbitrary set of input points (``testing_input``). It should be noted that, +# ``testing_trend`` should be set according to ``trend`` used during emulator +# training. + +testing_input = np.arange(0,10,0.1) +predictions = emulator.predict(testing_input) + +plt.plot(testing_input, predictions[:, 0], 'r-', label= "mean") +plt.scatter(x, y, s=15, c='k', label="training data", zorder=2) +plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI") +plt.xlabel('x') +plt.ylabel('emulator-predicted y') +plt.xlim(testing_input[0], testing_input[-1]) +_ = plt.legend() + +# %% md +# +# We can also draw any number of samples at ``testing_input``` using +# :py:meth:`.ScalarGaSP.sample()` method. + +nsamples = 5 +samples = emulator.sample(testing_input, nsamples=nsamples) + +for i in range(nsamples): + plt.plot(testing_input, samples[:,i], '--', label=f'sample{i+1}') + +plt.scatter(x, y, s=15, c='k', label="training data", zorder=2) +plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI") +plt.xlabel('x') +plt.ylabel('emulator-predicted y') +plt.xlim(testing_input[0], testing_input[-1]) +_ = plt.legend() \ No newline at end of file diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 2616ef0..9f07b14 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -3,4 +3,5 @@ Example Gallery .. toctree:: - Emulator - RobustGaSP <../auto_examples/emulator/robustgasp/index> \ No newline at end of file + Emulator - RobustGaSP <../auto_examples/emulator/robustgasp/index> + \ No newline at end of file -- GitLab