diff --git a/tests/AnalyticalFunction/Outputs_PostProcessing_calib/totalsobol_Z.csv b/tests/AnalyticalFunction/Outputs_PostProcessing_calib/totalsobol_Z.csv
index 51e13cf5fabf8c1d758cae33d1e7fe9e1deabe93..971b3917d08f2e8ad40102799b8ccaaee1b0ee7f 100644
--- a/tests/AnalyticalFunction/Outputs_PostProcessing_calib/totalsobol_Z.csv
+++ b/tests/AnalyticalFunction/Outputs_PostProcessing_calib/totalsobol_Z.csv
@@ -1,11 +1,11 @@
$X_1$,$X_2$
-9.993116948359226370e-01,3.945526636751626182e-01
-9.961447327121788486e-01,3.667186570442192428e-01
-9.947646363203489495e-01,3.647464239325610147e-01
-9.957567531711847275e-01,3.731759386892447572e-01
-9.938740724688709394e-01,3.496238556702422340e-01
-9.934295067779727040e-01,3.489816158549155500e-01
-1.000000000000000000e+00,2.733114696357644147e-01
-9.929163528111150105e-01,3.544359760567539208e-01
-9.927050241098982486e-01,3.584158152985525980e-01
-9.915077464334921542e-01,3.702678640580663316e-01
+9.366935781467614630e-01,1.296011206606640376e-01
+9.366778238387639677e-01,1.296333541870885453e-01
+9.366737110183974702e-01,1.296417560822205883e-01
+9.366715749203039509e-01,1.296461022957972054e-01
+9.366705838329740796e-01,1.296480821342860246e-01
+9.366707796164559507e-01,1.296475551494160605e-01
+9.366843916343090548e-01,1.296184596592687976e-01
+9.366660378655075059e-01,1.296580126549548295e-01
+9.366679119152595545e-01,1.296539068017917851e-01
+9.366689440746930329e-01,1.296517202892714382e-01
diff --git a/tests/AnalyticalFunction/AnalyticalFunction.py b/tests/AnalyticalFunction/analytical_function.py
similarity index 97%
rename from tests/AnalyticalFunction/AnalyticalFunction.py
rename to tests/AnalyticalFunction/analytical_function.py
index 0dbd4a97a4e000182e2f20afc8a535b79171e36d..7d1075a0195c7e3ab3db40465f02b840cebf28dc 100644
--- a/tests/AnalyticalFunction/AnalyticalFunction.py
+++ b/tests/AnalyticalFunction/analytical_function.py
@@ -11,7 +11,7 @@ import scipy.stats as st
import seaborn as sns
-def AnalyticalFunction(xx, t=None):
+def analytical_function(xx, t=None):
"""
Analytical Non-Gaussian Function
@@ -68,56 +68,56 @@ if __name__ == "__main__":
MCSize = 10000 #1000000
ndim = 10
sigma = 2
-
+
# -------------------------------------------------------------------------
# ----------------------- Synthetic data generation -----------------------
# -------------------------------------------------------------------------
t = np.arange(0, 10, 1.) / 9
-
+
MAP = np.zeros((1, ndim))
synthethicData = AnalyticalFunction(MAP, t=t)
-
+
# -------------------------------------------------------------------------
# ---------------------- Generate Prior distribution ----------------------
# -------------------------------------------------------------------------
-
+
xx = np.zeros((MCSize, ndim))
-
+
params = (-5,5)
-
+
for idxDim in range(ndim):
lower, upper = params
xx[:,idxDim] = stats.uniform(loc=lower, scale=upper-lower).rvs(size=MCSize)
-
+
# -------------------------------------------------------------------------
# ------------- BME and Kullback-Leibler Divergence -----------------------
# -------------------------------------------------------------------------
Outputs = AnalyticalFunction(xx, t=t)
-
+
cov_matrix = np.diag(np.repeat(sigma**2, synthethicData.shape[1]))
-
+
Likelihoods = st.multivariate_normal.pdf(Outputs[1:], mean=synthethicData[1], cov=cov_matrix)
-
+
sns.kdeplot(np.log(Likelihoods[Likelihoods>0]), shade=True, color="g", label='Ref. Likelihood')
-
+
normLikelihood = Likelihoods / np.nanmax(Likelihoods)
# Random numbers between 0 and 1
unif = np.random.rand(1, MCSize)[0]
-
+
# Reject the poorly performed prior
accepted = normLikelihood >= unif
-
+
# Prior-based estimation of BME
logBME = np.log(np.nanmean(Likelihoods))
print('\nThe Naive MC-Estimation of BME is %.5f.'%(logBME))
-
+
# Posterior-based expectation of likelihoods
postExpLikelihoods = np.mean(np.log(Likelihoods[accepted]))
-
+
# Calculate Kullback-Leibler Divergence
KLD = postExpLikelihoods - logBME
print("The Kullback-Leibler divergence estimation is %.5f."%KLD)
-
+
# -------------------------------------------------------------------------
# ----------------- Save the arrays as .npy files -------------------------
# -------------------------------------------------------------------------
diff --git a/tests/AnalyticalFunction/example_analytical_function.ipynb b/tests/AnalyticalFunction/example_analytical_function.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..e94174400f217ab2f8895c60ce37a2807b1a4864
--- /dev/null
+++ b/tests/AnalyticalFunction/example_analytical_function.ipynb
@@ -0,0 +1,890 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "13c8e941",
+ "metadata": {},
+ "source": [
+ "# Example: Surrogate Model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "de58892b",
+ "metadata": {},
+ "source": [
+ "In this example, we train a surrogate model via the `bayesvalidrox` package. We aim at training a Polynomial Chaos Expansion to a simple analytical function. \n",
+ "The PCE representation of the computational model $M$ provides the dependence of this model on the uncertain model's parameters $\\mathbf{\\theta}$ using projection onto an orthonormal polynomial basis. It could be also seen as a linear regression that includes linear combinations of a fixed set of nonlinear functions with respect to the input variables, known as polynomial basis function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2890d7ef",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:PCE_Trunc}\n",
+ " M(x,y,z,t, \\mathbf{\\theta}) \\approx \\sum_{\\mathbf{\\alpha} \\in \\mathcal{A} } c_{\\mathbf{\\alpha}} (x,y,z,t) \\Psi_{\\mathbf{\\alpha}}(\\mathbf{\\theta}) \\, .\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bccdaebc",
+ "metadata": {},
+ "source": [
+ "Here, $x,y,z,t$ are the spatial and temporal components of the quantity of interest, $\\mathbf{\\theta}$ is the vector of the $N$ uncertain parameters of model $M$, $c_{\\mathbf{\\alpha}}(x,y,z,t) \\in \\mathbb{R}$ are the corresponding expansion coefficients that are functions of space and time, and $\\Psi_{\\mathbf{\\alpha}}(\\mathbf{\\theta})$ represents multivariate polynomials orthogonal with respect to a multi-index $\\mathbf{\\alpha}$. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "50023bea",
+ "metadata": {},
+ "source": [
+ "The latter represents the combinatoric information how to enumerate all possible products of $N$ individual univariate basis functions with respect to the total degree of expansions less or equal to polynomial degree $d$:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "538c43c8",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:truncation}\n",
+ "\\begin{split}\n",
+ " \\mathcal{A}^{N, d} = \\{ \\alpha \\in \\mathbb{N}^{N} \\ : \\ |\\alpha|\\leq d\\} \\, , \\qquad\n",
+ " \\text{card} \\ \\mathcal{A}^{N, d} \\equiv P = \\binom{N+d}{d}.\n",
+ "\\end{split}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b22c0ccb",
+ "metadata": {},
+ "source": [
+ "The multivariate polynomials $\\Psi_{\\alpha}(\\mathbf{\\theta})$ are comprised of the tensor product of univariate polynomials"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "41eeb9d3",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:Psi}\n",
+ " \\Psi_{\\alpha}(\\mathbf{\\theta}_k) := \\prod_{i=1}^{N_k} \\psi_{\\alpha_i}^{(i)}(\\mathbf{\\theta}_{k,i}) \\, ,\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "860cbcbd",
+ "metadata": {},
+ "source": [
+ "where the univariate orthonormal polynomials $\\psi_{\\alpha_i}^{(i)}(\\mathbf{\\theta}_{i})$ must satisfy "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "459de9cc",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:univPsi}\n",
+ " \\langle \\psi_j^{(i)}(\\mathbf{\\theta}_{k,i}), \\psi_l^{(i)}(\\mathbf{\\theta}_{k,i}) \\rangle := \\int_{\\Theta_{k,i}} \\psi_j^{(i)}(\\mathbf{\\theta}_{k,i}) \\psi_l^{(i)}(\\mathbf{\\theta}_{k,i}) f_{\\Theta_{k,i}} (\\mathbf{\\theta}_{k,i})d \\mathbf{\\theta}_{k,i} = \\delta_{j l} \\, .\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "34403fc6",
+ "metadata": {},
+ "source": [
+ "Here, $i$ represents the input variable with respect to which the polynomials are orthogonal as well as the corresponding polynomial family, $j$ and $l$ are the corresponding polynomial degree, $f_{\\Theta_{i}}(\\mathbf{\\theta}_{i})$ is the $i$th-input marginal distribution and $\\delta_{j l}$ is the Kronecker delta.\n",
+ "We use an arbitrary polynomial chaos expansion (aPCE), introduced by [Oladyshkin & Nowak (2012)](https://www.sciencedirect.com/science/article/pii/S0951832012000853?casa_token=pbisUgY4niQAAAAA:8WsqMi1mCyfUIJ3GnFGdv6FXFA6a4g8MB75kjGGdEvocV64cd4E8LxcSh8_fwZTeI2ONlUalq_8), that can operate with probability measures that may be implicitly and incompletely defined via their statistical moments. Using aPCE, one can build the multivariate orthonormal polynomials even in the absence of the exact probability density function $f_{\\Theta}(\\theta)$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3bd3feba",
+ "metadata": {},
+ "source": [
+ "In this tutorial, we use an extension of aPCE as Bayesian sparse arbitrary polynomial chaos (BsaPCE) representation. This method computes the coefficients $c_\\alpha$ in a Bayesian setting via a so-called Bayesian sparse learning approach, introduced by [Tipping (2001)](https://www.jmlr.org/papers/volume1/tipping01a/tipping01a.pdf?ref=https://githubhelp.com)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6a291027",
+ "metadata": {},
+ "source": [
+ "The posterior distribution of the expansion coefficients, conditioned on the model responses $\\mathrm{\\mathbf{Y}}$ resulting from the training sets $\\mathbf{X}$, is given by the combination of a Gaussian likelihood and a Gaussian prior distribution over the unknown expansion coefficients $\\mathbf{c}$ according to Bayes' rule. Then, the posterior of the expansion coefficients given the model responses $\\mathrm{\\mathbf{Y}}$ and values of hyper-parameters $\\mathbf{\\alpha}$ and $\\beta$ describing the Gauss process, can take the following form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3127c49a",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:PCE_Posterior}\n",
+ " p(\\mathrm{\\mathbf{c}}|\\mathbf{Y},\\mathbf{\\alpha}, \\beta) = \\frac{p(\\mathrm{\\mathbf{Y}}|\\mathbf{X},\\mathbf{c}, \\beta) p(\\mathbf{c}|\\mathbf{\\alpha})}{p(\\mathrm{\\mathbf{Y}}| \\mathbf{X}, \\mathbf{\\alpha} , \\beta)},\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "46381c6d",
+ "metadata": {},
+ "source": [
+ "which is also Gaussian defined by $\\mathcal{N}( \\mathbf{c}| \\mathbf{\\mu}, \\mathbf{\\Sigma})$ with"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "515d47e6",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:PCE_Posterior_moments}\n",
+ " \\mathbf{\\mu} = \\beta \\mathbf{\\Sigma} \\mathbf{\\Psi}^{\\top} \\mathrm{\\mathbf{Y}} \\, , \\qquad\n",
+ " \\mathbf{\\Sigma} = \\left(\\mathbf{A}+ \\mathbf{\\Psi}^{\\top} \\beta \\mathbf{\\Psi} \\right)^{-1} \\, .\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ce4636d9",
+ "metadata": {},
+ "source": [
+ "Here, $\\mathbf{\\Psi}$ is the design matrix of size $E \\times N$ with elements $\\Psi_{ni}=\\psi_i(x_n)$, where $E$ represents the number of model evaluations using the training samples, and $\\mathbf{A}=\\mathrm{diag}(\\alpha_i)$. The values of $\\mathbf{\\alpha}$ and $\\beta$ can be determined via type-II maximum likelihood [Berger (2013)](https://books.google.com/books?hl=en&lr=&id=1CDaBwAAQBAJ&oi=fnd&pg=PA1&dq=Statistical+decision+theory+and+Bayesian+analysis.++berger+2013&ots=LMulrdTL3O&sig=xMTVRCVf5scWBLQi98BgUie5d-M)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c50a317d",
+ "metadata": {},
+ "source": [
+ "## Problem description: Analytical function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d13b9da3",
+ "metadata": {},
+ "source": [
+ "This test shows a surrogate-assisted Bayesian calibration of a time dependent non-linear analytical function of ten ($n=10$) uncertain parameters $\\omega=\\{\\omega_1, ..., \\omega_n\\}$, which reads as:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b715ea87",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\mathbf{y}(\\boldsymbol{\\omega}, t)=\\left(\\omega_{1}^{2}+\\omega_{2}-1\\right)^{2}+\\omega_{1}^{2}+0.1 \\omega_{1} \\exp \\left(\\omega_{2}\\right)-2 \\omega_{1} \\sqrt{0.5 t}+1+\\sum_{i=2}^{n} \\frac{\\omega_{i}^{3}}{i}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6625cdf6",
+ "metadata": {},
+ "source": [
+ "where the prior parameter distribution $p(\\omega)$ is considered to be independent and uniform with $\\omega_i \\sim \\mathcal{U} (-5, 5)$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "efe8f59c",
+ "metadata": {},
+ "source": [
+ "## Import necessary libraries"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "31b60d45",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import sys\n",
+ "import joblib\n",
+ "from IPython.display import IFrame"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6d44d5af",
+ "metadata": {},
+ "source": [
+ "## Define the model with PyLinkForwardModel"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bf3f4dfe",
+ "metadata": {},
+ "source": [
+ "We use the`PyLinkForwardModel`object for this purpose. Fistly, we are going to import the `bayesvalidrox` package and then, we instantiate the `PyLinkForwardModel` object."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "id": "25a3b65a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from bayesvalidrox import PyLinkForwardModel\n",
+ "Model = PyLinkForwardModel()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0452955f",
+ "metadata": {},
+ "source": [
+ "Next, we will pass the `link_type`, `name` and `py_file` variables to the `Model` object. Since the analytical function is implmented as a python function in a separate file, we only need to pass it's name (without `.py` extension) to the object variable `py_file`. Note that the function name in the python script should match that of the script. For models implemented as a separate python file, the `link_type` is `Function` to be given as a string. The `name` variable takes any user defined string."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "f26dbacd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Model.link_type = 'Function'\n",
+ "Model.py_file = 'analytical_function'\n",
+ "Model.name = 'AnalyticFunc'"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e12ef2ab",
+ "metadata": {},
+ "source": [
+ "The model output name is defined as follows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "c2778a83",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Model.Output.names = ['Z'] # As a list of strings"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d4106fb6",
+ "metadata": {},
+ "source": [
+ "**Bonus**: For this example, we have a Monte-Carlo reference solution for the first moements (mean and standard deviation) of the analytical function. The numpy (`*.npy`) files can be found in the `data\\` directory. We will discuss the first two moments with our estimate moments using the surrogate model. These values can be passed in a form of a dictionary to the object variable `mc_reference`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "9dc9e177",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Model.mc_reference = {}\n",
+ "Model.mc_reference['Time [s]'] = np.arange(0, 10, 1.) / 9\n",
+ "Model.mc_reference['mean'] = np.load(f\"data/mean_2.npy\")\n",
+ "Model.mc_reference['std'] = np.load(f\"data/std_2.npy\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c2c68eee",
+ "metadata": {},
+ "source": [
+ "## Define probablistic input model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "af7f7d8f",
+ "metadata": {},
+ "source": [
+ "Now, we define the distribution of the model inputs. `bayesvalidrox` accepts the definition in two ways: by defining the distribution directly or by passing available data. The latter is handy, when little information is available on the parameters or they do not follow any typical distribution. We will use the second option and read the input parameters form a numpy file in the `data/` directory."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "f1d2deb0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Import and instantiate the input object\n",
+ "from bayesvalidrox import Input\n",
+ "Inputs = Input()\n",
+ "\n",
+ "# Option I: Define distribution directy with their name, type and parameters \n",
+ "\n",
+ "# First parameter\n",
+ "# Inputs.addMarginals()\n",
+ "# Inputs.Marginals[0].name = '$X_1$'\n",
+ "# Inputs.Marginals[0].dist_type = 'unif'\n",
+ "# Inputs.Marginals[0].parameters = [-5, 5]\n",
+ "\n",
+ "# Second parameter\n",
+ "# Inputs.addMarginals()\n",
+ "# Inputs.Marginals[1].name = '$X_2$'\n",
+ "# Inputs.Marginals[1].dist_type = 'unif'\n",
+ "# Inputs.Marginals[1].parameters = [-5, 5]\n",
+ "\n",
+ "# ----------------------------------------------------------------------------\n",
+ "\n",
+ "# Option II: Pass available data for input parameters\n",
+ "inputParams = np.load('data/InputParameters_2.npy')\n",
+ "\n",
+ "# First parameter\n",
+ "Inputs.addMarginals()\n",
+ "Inputs.Marginals[0].name = '$X_1$'\n",
+ "Inputs.Marginals[0].input_data = inputParams[:, 0]\n",
+ "\n",
+ "# Second parameter\n",
+ "Inputs.addMarginals()\n",
+ "Inputs.Marginals[1].name = '$X_2$'\n",
+ "Inputs.Marginals[1].input_data = inputParams[:, 1]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "10ff4b23",
+ "metadata": {},
+ "source": [
+ "## Define surrogate (meta) model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d0e8a911",
+ "metadata": {},
+ "source": [
+ "In this example, we use a Polynomial Chaos Expansion (PCE) as our meta model. Like before, we need to import the `MetaModel` object from `bayesvalidrox` package and instantiate a meta-model object. This object, however, accepts the input object (`Input`) as an argument."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "id": "7256d8f0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from bayesvalidrox import MetaModel\n",
+ "MetaModelOpts = MetaModel(Inputs)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aa7756f5",
+ "metadata": {},
+ "source": [
+ "Let us define now the meta model type, the regression type, degree of the polynomials and the trunction norm `q_norm` which lies between 0 and 1. This parameter defines hyperbolic truncation scheme."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "id": "39a2ece3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Select your metamodel method\n",
+ "MetaModelOpts.meta_model_type = 'PCE'\n",
+ "\n",
+ "# Select the regression method for calculation of the PCE coefficients:\n",
+ "# 1)OLS: Ordinary Least Square 2)BRR: Bayesian Ridge Regression\n",
+ "# 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression\n",
+ "# 5)FastARD: Fast Bayesian ARD Regression\n",
+ "# 6)VBL: Variational Bayesian Learning\n",
+ "# 7)EBL: Emperical Bayesian Learning\n",
+ "MetaModelOpts.pce_reg_method = 'FastARD'\n",
+ "\n",
+ "\n",
+ "# Specify the polynomial degree to be compared by the adaptive algorithm:\n",
+ "# The degree with the lowest Leave-One-Out cross-validation (LOO)\n",
+ "# error (or the highest score=1-LOO)estimator is chosen as the final\n",
+ "# metamodel. pce_deg accepts degree as a scalar or a range.\n",
+ "MetaModelOpts.pce_deg = np.arange(9)\n",
+ "\n",
+ "# Hyperbolic truncation scheme 0<q<1 (default=1)\n",
+ "MetaModelOpts.pce_q_norm = 0.75"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a3048f38",
+ "metadata": {},
+ "source": [
+ "After defining the metamodel type, we need to define the so-called experimental design (ExpDesign). ExpDesign basically provides instruction on how to samplie the input parameter space."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "id": "f8e1bf2d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ------------------------------------------------\n",
+ "# ------ Experimental Design Configuration -------\n",
+ "# ------------------------------------------------\n",
+ "MetaModelOpts.add_ExpDesign()\n",
+ "\n",
+ "# One-shot (normal) or Sequential Adaptive (sequential) Design\n",
+ "MetaModelOpts.ExpDesign.Method = 'normal'\n",
+ "MetaModelOpts.ExpDesign.n_init_samples = 100\n",
+ "\n",
+ "# Sampling methods\n",
+ "# 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) korobov\n",
+ "# 7) chebyshev(FT) 8) grid(FT) 9) nested_grid(FT) 10)user\n",
+ "MetaModelOpts.ExpDesign.sampling_method = 'latin_hypercube'"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8f7cd8dc",
+ "metadata": {},
+ "source": [
+ "Now, we can start training the surrogate (meta-) model by using the `create_metamodel` method and passing the model object as the only argument."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "id": "82fffca5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Computing orth. polynomial coeffs: 100%|##########| 2/2 [00:00<00:00, 3.42it/s]"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ " Now the forward model needs to be run!\n",
+ "\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "Running forward model : 100%|██████████| 100/100 [00:00<00:00, 6072.54it/s]\n",
+ "Fitting regression: 0%| | 0/1 [00:00<?, ?it/s]"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ ">>>> Training the PCE metamodel started. <<<<<<\n",
+ "\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Fitting regression: 100%|██████████| 1/1 [00:03<00:00, 3.76s/it]"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ ">>>> Training the PCE metamodel sucessfully completed. <<<<<<\n",
+ "\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Train the meta model\n",
+ "PCEModel = MetaModelOpts.create_metamodel(Model)\n",
+ "\n",
+ "# Save PCE models as pkl object for further deployment for example on a cloud\n",
+ "with open(f'PCEModel_{Model.name}.pkl', 'wb') as output:\n",
+ " joblib.dump(PCEModel, output, 2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4b2b6a1b",
+ "metadata": {},
+ "source": [
+ "## Post-processing"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "64ea5912",
+ "metadata": {},
+ "source": [
+ "As before, we need to import the `PostProcessing` module of `bayesvalidrox` and instantiate it. Bear in mind that it accepts the meta-model object `PCEModel` as the only argument."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "id": "d51e4195",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from bayesvalidrox import PostProcessing\n",
+ "PostPCE = PostProcessing(PCEModel)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "51aed43d",
+ "metadata": {},
+ "source": [
+ "### Moment comparison"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "69de7856",
+ "metadata": {},
+ "source": [
+ "Since the reference moments obtained from a Monte-Carlo simulation is available, we only need to call the `plotMoments` method from the PostProcessing object. This method generates a plot and stores it in `Outputs_PostProcessing_calib` directory."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "id": "96e85fb9",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ " <iframe\n",
+ " width=\"600\"\n",
+ " height=\"300\"\n",
+ " src=\"./Outputs_PostProcessing_calib/Mean_Std_PCE.pdf\"\n",
+ " frameborder=\"0\"\n",
+ " allowfullscreen\n",
+ " ></iframe>\n",
+ " "
+ ],
+ "text/plain": [
+ "<IPython.lib.display.IFrame at 0x7f88751bc070>"
+ ]
+ },
+ "execution_count": 34,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "text/plain": [
+ "<Figure size 1728x1152 with 0 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Compute the moments and compare with the Monte-Carlo reference\n",
+ "PostPCE.plotMoments()\n",
+ "# Show the pdf\n",
+ "IFrame(\"./Outputs_PostProcessing_calib/Mean_Std_PCE.pdf\", width=900, height=600)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a6cf5333",
+ "metadata": {},
+ "source": [
+ "### Validation of the metamodel"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9d506ae2",
+ "metadata": {},
+ "source": [
+ "Let us first visually compare the results of the metamodel and the original model, i.e. `Analyrical Function` for 3 randomly drawn samples for the prior parameter distribution. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "id": "cfd42c02",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Running forward model valid: 100%|██████████| 3/3 [00:00<00:00, 25627.11it/s]\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ " <iframe\n",
+ " width=\"600\"\n",
+ " height=\"300\"\n",
+ " src=\"./Outputs_PostProcessing_calib/Model_vs_PCEModel.pdf\"\n",
+ " frameborder=\"0\"\n",
+ " allowfullscreen\n",
+ " ></iframe>\n",
+ " "
+ ],
+ "text/plain": [
+ "<IPython.lib.display.IFrame at 0x7f883b4a4490>"
+ ]
+ },
+ "execution_count": 37,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "text/plain": [
+ "<Figure size 1728x1152 with 0 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Plot to check validation visually.\n",
+ "PostPCE.validMetamodel(nValidSamples=3)\n",
+ "# Show the pdf\n",
+ "IFrame(\"./Outputs_PostProcessing_calib/Model_vs_PCEModel.pdf\", width=900, height=600)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3b39e05f",
+ "metadata": {},
+ "source": [
+ "Another way to check the accuracy of the meta model is to use the `accuracyCheckMetaModel` method to show the Root Mean Square Error and the validation error. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "id": "e2f83de2",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Running forward model validSet: 100%|██████████| 200/200 [00:00<00:00, 6802.31it/s]"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ ">>>>> Errors of Z <<<<<\n",
+ "\n",
+ "Index | RMSE | Validation Error\n",
+ "-----------------------------------\n",
+ "1 | 4.480e-01 | 4.290e-08\n",
+ "2 | 4.481e-01 | 4.289e-08\n",
+ "3 | 4.482e-01 | 4.288e-08\n",
+ "4 | 4.483e-01 | 4.289e-08\n",
+ "5 | 4.484e-01 | 4.292e-08\n",
+ "6 | 4.490e-01 | 4.301e-08\n",
+ "7 | 4.583e-01 | 4.481e-08\n",
+ "8 | 4.463e-01 | 4.249e-08\n",
+ "9 | 4.472e-01 | 4.264e-08\n",
+ "10 | 4.474e-01 | 4.268e-08\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Compute and print RMSE error\n",
+ "PostPCE.accuracyCheckMetaModel(nSamples=200)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ef686ae7",
+ "metadata": {},
+ "source": [
+ "### Global sensitivity analysis with Sobol indices"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1ffecd15",
+ "metadata": {},
+ "source": [
+ "Here, we analyze how the variability of the model response quantity (`Z`) is affected by the variability of each input variable or combinations thereof. Here, we use the so-called Sobol indices ([Sobol original paper](https://mae.ufl.edu/haftka/eoed/protected/Sobol%20Original%20Paper.pdf)), derived from a variance decomposition of model outputs in terms of contributions of each input parameter or combinations thereof. \n",
+ "Using Sobol decomposition, one can describe the total variance of the model in terms of the sum of the summands' variances. Once the PC representation of the model is available, the expansion coefficients are simply gathered according to the dependency of each basis polynomial, square-summed and normalized"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3e95b594",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:pce-sobol-1st}\n",
+ "\\begin{array}{l}\n",
+ "S_{i_{1}, \\ldots, i_{s}}=\\frac{\\sum\\limits_{j=1}^{M} \\chi_{j} c_{j}^{2}}{\\sum\\limits_{j=1}^{M} c_{j}^{2}} \\, ,\\qquad\n",
+ "\\chi_{j}=\\left\\{\\begin{array}{ll}\n",
+ "1, & \\text { if } \\alpha_{j}^{k}>0, \\quad \\forall j \\in\\left(i_{1}, \\ldots, i_{s}\\right) \\\\[0.5em]\n",
+ "0, & \\text { if } \\alpha_{j}^{k}=0, \\quad \\exists j \\in\\left(i_{1}, \\ldots, i_{s}\\right)\n",
+ "\\end{array}\\right\\} \\, .\n",
+ "\\end{array}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a255fc43",
+ "metadata": {},
+ "source": [
+ "Here, $S_{i_{1}, \\ldots, i_{s}}$ is the Sobol index that indicates what fraction of total variance of the response quantity can be traced back to the joint contributions of the parameters $\\theta_{i_{1}}, \\ldots, \\theta_{i_{s}}.$ The index selection operator $\\chi_{j}$ indicates where the chosen parameters $\\theta$ numbered as $i_{1}, \\ldots, i_{s}$ (i.e., $\\left.\\theta_{i_{1}}, \\ldots, \\theta_{i_{s}}\\right)$ have concurrent contributions to the variance within the overall expansion. Simply put, it selects all polynomial terms with the specified combination $i_{1}, \\ldots, i_{s}$ of model parameters."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "025e4a83",
+ "metadata": {},
+ "source": [
+ "A complementing measure for sensitivity analysis is the Sobol Total Index. It expresses the total contribution to the variance of model output due to the uncertainty of an individual parameter $\\theta_j$ in all cross-combinations with other parameters"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9da5d910",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\label{eq:pce-sobol-total}\n",
+ "S_{j}^{T}=\\sum_{\\left\\{i_{1}, \\ldots, i_{s}\\right\\} \\supset j} S_{i_{1}, \\ldots, i_{s}},\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "59413b8c",
+ "metadata": {},
+ "source": [
+ "where $S_{j}^{T}$ is simply a summation of all Sobol indices in which the variable $\\theta_j$ appears as univariate as well as joint influences.\n",
+ "The Total Sobol indices sum to one, if input variables are independent. When dealing with correlated variables, however, this is not the case."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "962bf769",
+ "metadata": {},
+ "source": [
+ "To perform the sensitivity analysis with `bayesvalidrox` package, we need to call the `sobolIndicesPCE` of the `PostProcessing` object. This returns two dictionaries containing the single sobol indices and the total ones. Moreover, it plots the Total Sobol Indices and stores the plots in `pdf` format in `Outputs_PostProcessing_calib` directory."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 42,
+ "id": "135173d6",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ " <iframe\n",
+ " width=\"900\"\n",
+ " height=\"600\"\n",
+ " src=\"./Outputs_PostProcessing_calib/Sobol_indices.pdf\"\n",
+ " frameborder=\"0\"\n",
+ " allowfullscreen\n",
+ " ></iframe>\n",
+ " "
+ ],
+ "text/plain": [
+ "<IPython.lib.display.IFrame at 0x7f883b4cab80>"
+ ]
+ },
+ "execution_count": 42,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "text/plain": [
+ "<Figure size 1728x1152 with 0 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Plot the sobol indices\n",
+ "sobol_cell, total_sobol = PostPCE.sobolIndicesPCE()\n",
+ "# Show the pdf\n",
+ "IFrame(\"./Outputs_PostProcessing_calib/Sobol_indices.pdf\", width=900, height=600)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}