diff --git a/examples/borehole/borehole.py b/examples/borehole/borehole.py
new file mode 100644
index 0000000000000000000000000000000000000000..43dd71f21904a7ca7bec4e45a48fd3a37cc33879
--- /dev/null
+++ b/examples/borehole/borehole.py
@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 19 08:56:21 2018
+
+@author: farid
+"""
+import numpy as np
+
+
+def borehole(xx, *args):
+ """
+ BOREHOLE FUNCTION
+
+ Authors: Sonja Surjanovic, Simon Fraser University
+ Derek Bingham, Simon Fraser University
+ Questions/Comments: Please email Derek Bingham at dbingham@stat.sfu.ca.
+
+ Copyright 2013. Derek Bingham, Simon Fraser University.
+
+ THERE IS NO WARRANTY, EXPRESS OR IMPLIED. WE DO NOT ASSUME ANY LIABILITY
+ FOR THE USE OF THIS SOFTWARE. If software is modified to produce
+ derivative works, such modified software should be clearly marked.
+ Additionally, this program is free software; you can redistribute it
+ and/or modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; version 2.0 of the License.
+ Accordingly, this program is distributed in the hope that it will be
+ useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ For function details and reference information, see:
+ https://www.sfu.ca/~ssurjano/ishigami.html
+
+ Parameters
+ ----------
+ xx : array of shape (n_samples, n_params)
+ Input parameter sets. [rw, L, kw, Tu, Tl, Hu, Hl, r]
+
+ Returns
+ -------
+ output: dict
+ Water flow rate.
+
+ """
+ rw = xx[:, 0]
+ L = xx[:, 1]
+ Kw = xx[:, 2]
+ Tu = xx[:, 3]
+ Tl = xx[:, 4]
+ Hu = xx[:, 5]
+ Hl = xx[:, 6]
+ r = xx[:, 7]
+
+ frac1 = 2 * np.pi * Tu * (Hu-Hl)
+
+ frac2a = 2*L*Tu / (np.log(r/rw)*(rw**2)*Kw)
+ frac2b = Tu / Tl
+ frac2 = np.log(r/rw) * (1 + frac2a + frac2b)
+
+ y = frac1 / frac2
+
+ # Prepare output dict using standard bayesvalidrox format
+ output = {'x_values': np.zeros(1), 'flow rate [m$^3$/yr]': y}
+
+ return output
diff --git a/examples/borehole/data/valid_outputs.npy b/examples/borehole/data/valid_outputs.npy
new file mode 100644
index 0000000000000000000000000000000000000000..076320cebb087989776523cbfd6e1dd122c84b63
Binary files /dev/null and b/examples/borehole/data/valid_outputs.npy differ
diff --git a/examples/borehole/data/valid_samples.npy b/examples/borehole/data/valid_samples.npy
new file mode 100644
index 0000000000000000000000000000000000000000..27564f822ee863e5d75d7b8c71fcd8a763d5bc63
Binary files /dev/null and b/examples/borehole/data/valid_samples.npy differ
diff --git a/examples/borehole/test_borehole.py b/examples/borehole/test_borehole.py
new file mode 100644
index 0000000000000000000000000000000000000000..3097b78bc8fd221d18f254b2426b9fdff00987e8
--- /dev/null
+++ b/examples/borehole/test_borehole.py
@@ -0,0 +1,240 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This test deals with the surrogate modeling of a Ishigami function.
+
+You will see how to:
+ Check the quality of your regression model
+ Perform sensitivity analysis via Sobol Indices
+
+Author: Farid Mohammadi, M.Sc.
+E-Mail: farid.mohammadi@iws.uni-stuttgart.de
+Department of Hydromechanics and Modelling of Hydrosystems (LH2)
+Institute for Modelling Hydraulic and Environmental Systems (IWS), University
+of Stuttgart, www.iws.uni-stuttgart.de/lh2/
+Pfaffenwaldring 61
+70569 Stuttgart
+
+Created on Wed Jul 10 2019
+
+"""
+
+import numpy as np
+import joblib
+
+# import bayesvalidrox
+# Add BayesValidRox path
+import sys
+sys.path.append("../../src/bayesvalidrox/")
+
+from pylink.pylink import PyLinkForwardModel
+from surrogate_models.inputs import Input
+from surrogate_models.surrogate_models import MetaModel
+from post_processing.post_processing import PostProcessing
+from bayes_inference.bayes_inference import BayesInference
+from bayes_inference.discrepancy import Discrepancy
+import matplotlib
+matplotlib.use('agg')
+
+if __name__ == "__main__":
+
+ # =====================================================
+ # ============= COMPUTATIONAL MODEL ================
+ # =====================================================
+ Model = PyLinkForwardModel()
+
+ # Define model options
+ Model.link_type = 'Function'
+ Model.py_file = 'borehole'
+ Model.name = 'borehole'
+
+ Model.Output.names = ['flow rate [m$^3$/yr]']
+
+ # =====================================================
+ # ========= PROBABILISTIC INPUT MODEL ==============
+ # =====================================================
+ Inputs = Input()
+
+ # borehole radius
+ Inputs.add_marginals()
+ Inputs.Marginals[0].name = '$r_w$'
+ Inputs.Marginals[0].dist_type = 'normal'
+ Inputs.Marginals[0].parameters = [0.10, 0.0161812]
+
+ # borehole length
+ Inputs.add_marginals()
+ Inputs.Marginals[1].name = '$L$'
+ Inputs.Marginals[1].dist_type = 'unif'
+ Inputs.Marginals[1].parameters = [1120, 1680]
+
+ # borehole hydraulic conductivity
+ Inputs.add_marginals()
+ Inputs.Marginals[2].name = '$K_w$'
+ Inputs.Marginals[2].dist_type = 'unif'
+ Inputs.Marginals[2].parameters = [9855, 12045]
+
+ # transmissivity of upper aquifer
+ Inputs.add_marginals()
+ Inputs.Marginals[3].name = '$T_u$'
+ Inputs.Marginals[3].dist_type = 'unif'
+ Inputs.Marginals[3].parameters = [63070, 115600]
+
+ # transmissivity of lower aquifer
+ Inputs.add_marginals()
+ Inputs.Marginals[4].name = '$T_l$'
+ Inputs.Marginals[4].dist_type = 'unif'
+ Inputs.Marginals[4].parameters = [63.1, 116]
+
+ # potentiometric head of upper aquifer
+ Inputs.add_marginals()
+ Inputs.Marginals[5].name = '$H_u$'
+ Inputs.Marginals[5].dist_type = 'unif'
+ Inputs.Marginals[5].parameters = [990, 1110]
+
+ # potentiometric head of lower aquifer
+ Inputs.add_marginals()
+ Inputs.Marginals[6].name = '$H_l$'
+ Inputs.Marginals[6].dist_type = 'unif'
+ Inputs.Marginals[6].parameters = [700, 820]
+
+ # radius of influence
+ Inputs.add_marginals()
+ Inputs.Marginals[7].name = '$r$'
+ Inputs.Marginals[7].dist_type = 'lognorm'
+ Inputs.Marginals[7].parameters = [7.71, 1.0056]
+
+ # =====================================================
+ # ====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
+ # =====================================================
+ MetaModelOpts = MetaModel(Inputs, Model)
+
+ # Select your metamodel method
+ # 1) PCE (Polynomial Chaos Expansion) 2) aPCE (arbitrary PCE)
+ # 3) GPE (Gaussian Process Emulator)
+ MetaModelOpts.meta_model_type = 'aPCE'
+
+ # ------------------------------------------------
+ # ------------- PCE Specification ----------------
+ # ------------------------------------------------
+ # Select the sparse least-square minimization method for
+ # the PCE coefficients calculation:
+ # 1)OLS: Ordinary Least Square 2)BRR: Bayesian Ridge Regression
+ # 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression
+ # 5)FastARD: Fast Bayesian ARD Regression
+ # 6)BCS: Bayesian Compressive Sensing
+ # 7)OMP: Orthogonal Matching Pursuit
+ # 8)VBL: Variational Bayesian Learning
+ # 9)EBL: Emperical Bayesian Learning
+ MetaModelOpts.pce_reg_method = 'OMP'
+
+ # Bootstraping
+ # 1) normal 2) fast
+ MetaModelOpts.bootstrap_method = 'fast'
+ MetaModelOpts.n_bootstrap_itrs = 1#00
+
+ # Specify the max degree to be compared by the adaptive algorithm:
+ # The degree with the lowest Leave-One-Out cross-validation (LOO)
+ # error (or the highest score=1-LOO)estimator is chosen as the final
+ # metamodel. pce_deg accepts degree as a scalar or a range.
+ MetaModelOpts.pce_deg = 5
+
+ # q-quasi-norm 0<q<1 (default=1)
+ MetaModelOpts.pce_q_norm = 1.0
+
+ # Print summary of the regression results
+ # MetaModelOpts.verbose = True
+
+ # ------------------------------------------------
+ # ------ Experimental Design Configuration -------
+ # ------------------------------------------------
+ MetaModelOpts.add_ExpDesign()
+
+ # One-shot (normal) or Sequential Adaptive (sequential) Design
+ MetaModelOpts.ExpDesign.method = 'sequential'
+ MetaModelOpts.ExpDesign.n_init_samples = 50
+
+ # Sampling methods
+ # 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) korobov
+ # 7) chebyshev(FT) 8) grid(FT) 9) nested_grid(FT) 10)user
+ MetaModelOpts.ExpDesign.sampling_method = 'latin_hypercube'
+
+ # Provide the experimental design object with a hdf5 file
+ # MetaModelOpts.ExpDesign.hdf5_file = 'ExpDesign_borehole.hdf5'
+
+ # Sequential experimental design (needed only for sequential ExpDesign)
+ MetaModelOpts.ExpDesign.n_new_samples = 1
+ MetaModelOpts.ExpDesign.n_max_samples = 300 # 150
+ MetaModelOpts.ExpDesign.mod_LOO_threshold = 1e-16
+
+ # ------------------------------------------------
+ # ------- Sequential Design configuration --------
+ # ------------------------------------------------
+ # 1) None 2) 'equal' 3)'epsilon-decreasing' 4) 'adaptive'
+ MetaModelOpts.ExpDesign.tradeoff_scheme = None
+ # MetaModelOpts.ExpDesign.n_replication = 50
+ # -------- Exploration ------
+ # 1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'dual annealing'
+ MetaModelOpts.ExpDesign.explore_method = 'latin_hypercube'
+
+ # Use when 'dual annealing' chosen
+ MetaModelOpts.ExpDesign.max_func_itr = 200
+
+ # Use when 'Voronoi' or 'random' or 'latin_hypercube' chosen
+ MetaModelOpts.ExpDesign.n_canddidate = 5000
+ MetaModelOpts.ExpDesign.n_cand_groups = 4
+
+ # -------- Exploitation ------
+ # 1)'BayesOptDesign' 2)'VarOptDesign' 3)'alphabetic' 4)'Space-filling'
+ MetaModelOpts.ExpDesign.exploit_method = 'Space-filling'
+
+ # BayesOptDesign -> when data is available
+ # 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
+ # 3)APP (A-Posterior-percision)
+ # MetaModelOpts.ExpDesign.util_func = 'DKL'
+
+ # VarBasedOptDesign -> when data is not available
+ # Only with Vornoi >>> 1)Entropy 2)EIGF, 3)ALM, 4)LOOCV
+ MetaModelOpts.ExpDesign.util_func = 'ALM'
+
+ # alphabetic
+ # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
+ # 3)K-Opt (K-Optimality)
+ # MetaModelOpts.ExpDesign.util_func = 'D-Opt'
+
+ MetaModelOpts.valid_samples = np.load("data/valid_samples.npy")
+ MetaModelOpts.valid_model_runs = {
+ 'flow rate [m$^3$/yr]': np.load("data/valid_outputs.npy")
+ }
+ # >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ # Adaptive sparse arbitrary polynomial chaos expansion
+ # PCEModel = MetaModelOpts.create_metamodel()
+ from surrogate_models.meta_model_engine import MetaModelEngine
+ meta_model_engine = MetaModelEngine(MetaModelOpts)
+ meta_model_engine.run()
+ PCEModel = meta_model_engine.MetaModel
+
+ # Save PCE models
+ with open(f'PCEModel_{Model.name}.pkl', 'wb') as output:
+ joblib.dump(PCEModel, output, 2)
+
+ # =====================================================
+ # ========= POST PROCESSING OF METAMODELS ===========
+ # =====================================================
+ PostPCE = PostProcessing(PCEModel)
+
+ # Plot to check validation visually.
+ PostPCE.valid_metamodel(n_samples=200)
+
+ # Check the quality of your regression model
+ PostPCE.check_reg_quality()
+
+ # PostPCE.eval_PCEmodel_3D()
+ # Compute and print RMSE error
+ PostPCE.check_accuracy(n_samples=3000)
+
+ # Plot the evolution of the KLD,BME, and Modified LOOCV error
+ if MetaModelOpts.ExpDesign.method == 'sequential':
+ PostPCE.plot_seq_design_diagnostics()
+
+ # Plot the sobol indices
+ total_sobol = PostPCE.sobol_indices(plot_type='bar')