[examples] add borehole example.

examples/borehole/borehole.py

+#!/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

examples/borehole/test_borehole.py

+#!/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')