From 671fad7e02364ed9198c9b3bf686de63a73aa6cd Mon Sep 17 00:00:00 2001
From: Farid Mohammadi <farid.mohammadi@iws.uni-stuttgart.de>
Date: Wed, 14 Sep 2022 16:08:40 +0200
Subject: [PATCH] [examples] add OHagan example.
---
examples/OHagan-function/OHagan.py | 58 +++++
examples/OHagan-function/data/M.csv | 15 ++
.../data/sparse_solver_comparison.py | 240 ++++++++++++++++++
.../OHagan-function/data/valid_outputs.npy | Bin 0 -> 928 bytes
.../OHagan-function/data/valid_samples.npy | Bin 0 -> 12128 bytes
examples/OHagan-function/test_OHagan.py | 197 ++++++++++++++
6 files changed, 510 insertions(+)
create mode 100644 examples/OHagan-function/OHagan.py
create mode 100644 examples/OHagan-function/data/M.csv
create mode 100644 examples/OHagan-function/data/sparse_solver_comparison.py
create mode 100644 examples/OHagan-function/data/valid_outputs.npy
create mode 100644 examples/OHagan-function/data/valid_samples.npy
create mode 100644 examples/OHagan-function/test_OHagan.py
diff --git a/examples/OHagan-function/OHagan.py b/examples/OHagan-function/OHagan.py
new file mode 100644
index 000000000..71a58e092
--- /dev/null
+++ b/examples/OHagan-function/OHagan.py
@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 19 08:56:21 2018
+
+@author: farid
+"""
+import numpy as np
+
+
+def OHagan(xx, *args):
+ """
+ Oakley & O'Hagan (2004) Function
+
+ This function's a-coefficients are chosen so that 5 of the input variables
+ contribute significantly to the output variance, 5 have a much smaller
+ effect, and the remaining 5 have almost no effect on the output variance.
+
+ O'Hagan, 2004, Probabilistic sensitivity analysis of complex models: a
+ Bayesian approach J. R. Statist. Soc. B (2004) 66, Part 3, pp. 751-769.
+
+ Parameters
+ ----------
+ xx : array of shape (n_samples, n_params)
+ Input parameter sets.
+
+ Returns
+ -------
+ output: dict
+ Ourput value.
+
+ """
+ n_samples, n_params = xx.shape
+
+ # Read M
+ M = np.loadtxt('data/M.csv', delimiter=',')
+
+ a_1 = np.array([0.0118, 0.0456, 0.2297, 0.0393, 0.1177, 0.3865, 0.3897,
+ 0.6061, 0.6159, 0.4005, 1.0741, 1.1474, 0.7880, 1.1242,
+ 1.1982])
+ a_2 = np.array([0.4341, 0.0887, 0.0512, 0.3233, 0.1489, 1.0360, 0.9892,
+ 0.9672, 0.8977, 0.8083, 1.8426, 2.4712, 2.3946, 2.0045,
+ 2.2621])
+ a_3 = np.array([0.1044, 0.2057, 0.0774, 0.2730, 0.1253, 0.7526, 0.8570,
+ 1.0331, 0.8388, 0.7970, 2.2145, 2.0382, 2.4004, 2.0541,
+ 1.9845])
+
+ # Compute response
+ xx_M = np.dot(xx, M)
+ last_term = np.sum(xx_M * xx, axis=1)
+
+ y = np.dot(xx, a_1) + np.dot(np.sin(xx), a_2) + np.dot(np.cos(xx), a_3)
+ y += last_term
+
+ # Prepare output dict using standard bayesvalidrox format
+ output = {'x_values': np.zeros(1), 'Z': y.T}
+
+ return output
diff --git a/examples/OHagan-function/data/M.csv b/examples/OHagan-function/data/M.csv
new file mode 100644
index 000000000..98ec96817
--- /dev/null
+++ b/examples/OHagan-function/data/M.csv
@@ -0,0 +1,15 @@
+-0.022482886,-0.18501666,0.13418263,0.36867264,0.17172785,0.13651143,-0.44034404,-0.081422854,0.71321025,-0.44361072,0.50383394,-0.024101458,-0.045939684,0.21666181,0.055887417
+0.2565963,0.053792287,0.25800381,0.23795905,-0.59125756,-0.081627077,-0.28749073,0.41581639,0.49752241,0.083893165,-0.11056683,0.033222351,-0.13979497,-0.031020556,-0.22318721
+-0.055999811,0.19542252,0.095529005,-0.2862653,-0.14441303,0.22369356,0.14527412,0.28998481,0.2310501,-0.31929879,-0.29039128,-0.20956898,0.43139047,0.024429152,0.044904409
+0.66448103,0.43069872,0.29924645,-0.16202441,-0.31479544,-0.39026802,0.17679822,0.057952663,0.17230342,0.13466011,-0.3527524,0.25146896,-0.018810529,0.36482392,-0.32504618
+-0.121278,0.12463327,0.10656519,0.046562296,-0.21678617,0.19492172,-0.065521126,0.024404669,-0.09682886,0.19366196,0.33354757,0.31295994,-0.083615456,-0.25342082,0.37325717
+-0.2837623,-0.32820154,-0.10496068,-0.22073452,-0.13708154,-0.14426375,-0.11503319,0.22424151,-0.030395022,-0.51505615,0.017254978,0.038957118,0.36069184,0.30902452,0.050030193
+-0.077875893,0.003745656,0.88685604,-0.26590028,-0.079325357,-0.042734919,-0.18653782,-0.35604718,-0.17497421,0.088699956,0.40025886,-0.055979693,0.13724479,0.21485613,-0.011265799
+-0.09229473,0.59209563,0.031338285,-0.033080861,-0.24308858,-0.099798547,0.034460195,0.095119813,-0.3380162,0.0063860024,-0.61207299,0.081325416,0.88683114,0.14254905,0.14776204
+-0.13189434,0.52878496,0.12652391,0.045113625,0.58373514,0.37291503,0.11395325,-0.29479222,-0.57014085,0.46291592,-0.094050179,0.13959097,-0.38607402,-0.4489706,-0.14602419
+0.058107658,-0.32289338,0.093139162,0.072427234,-0.56919401,0.52554237,0.23656926,-0.011782016,0.071820601,0.078277291,-0.13355752,0.22722721,0.14369455,-0.45198935,-0.55574794
+0.66145875,0.34633299,0.14098019,0.51882591,-0.28019898,-0.1603226,-0.068413337,-0.20428242,0.069672173,0.23112577,-0.044368579,-0.16455425,0.21620977,0.0042702105,-0.087399014
+0.31599556,-0.027551859,0.13434254,0.13497371,0.05400568,-0.17374789,0.17525393,0.060258929,-0.17914162,-0.31056619,-0.25358691,0.025847535,-0.43006001,-0.62266361,-0.033996882
+-0.29038151,0.03410127,0.034903413,-0.12121764,0.026030714,-0.33546274,-0.41424111,0.05324838,-0.27099455,-0.026251302,0.41024137,0.26636349,0.15582891,-0.18666254,0.019895831
+-0.24388652,-0.44098852,0.012618825,0.24945112,0.071101888,0.24623792,0.17484502,0.0085286769,0.2514707,-0.14659862,-0.08462515,0.36931333,-0.29955293,0.1104436,-0.75690139
+0.041494323,-0.25980564,0.46402128,-0.36112127,-0.94980789,-0.16504063,0.0030943325,0.052792942,0.22523648,0.38390366,0.45562427,-0.18631744,0.0082333995,0.16670803,0.16045688
diff --git a/examples/OHagan-function/data/sparse_solver_comparison.py b/examples/OHagan-function/data/sparse_solver_comparison.py
new file mode 100644
index 000000000..5245bb89d
--- /dev/null
+++ b/examples/OHagan-function/data/sparse_solver_comparison.py
@@ -0,0 +1,240 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Sep 10 09:44:05 2022
+
+@author: farid
+"""
+
+import numpy as np
+import joblib
+import os
+import scipy.stats as st
+
+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')
+from matplotlib.backends.backend_pdf import PdfPages
+import matplotlib.ticker as ticker
+from matplotlib.offsetbox import AnchoredText
+from matplotlib.patches import Patch
+import matplotlib.pyplot as plt
+# Load the mplstyle
+plt.style.use(os.path.join(
+ os.path.split(__file__)[0],
+ '../../../src/bayesvalidrox/', 'bayesvalidrox.mplstyle'))
+
+
+def plot_seq_design_diagnostics(meta_model, util_funcs,
+ ref_BME_KLD=None, save_fig=True):
+ """
+ Plots the Bayesian Model Evidence (BME) and Kullback-Leibler divergence
+ (KLD) for the sequential design.
+
+ Parameters
+ ----------
+ ref_BME_KLD : array, optional
+ Reference BME and KLD . The default is `None`.
+ save_fig : bool, optional
+ Whether to save the figures. The default is `True`.
+
+ Returns
+ -------
+ None.
+
+ """
+ n_init_samples = meta_model.ExpDesign.n_init_samples
+ n_total_samples = meta_model.ExpDesign.X.shape[0]
+
+ if save_fig:
+ newpath = f'{path}/boxplot_{model_name}/'
+ if not os.path.exists(newpath):
+ os.makedirs(newpath)
+
+ plotList = ['Modified LOO error', 'Validation error', 'KLD', 'BME',
+ 'RMSEMean', 'RMSEStd', 'Hellinger distance']
+ seqList = [meta_model.SeqModifiedLOO, meta_model.seqValidError,
+ meta_model.SeqKLD, meta_model.SeqBME, meta_model.seqRMSEMean,
+ meta_model.seqRMSEStd, meta_model.SeqDistHellinger]
+
+ markers = ('x', 'o', 'd', '*', '+')
+ colors = ('k', 'darkgreen', 'b', 'navy', 'darkred')
+
+ # Plot the evolution of the diagnostic criteria of the
+ # Sequential Experimental Design.
+ for plotidx, plot in enumerate(plotList):
+ fig, ax = plt.subplots(figsize=(27, 15))
+ seq_dict = seqList[plotidx]
+ name_util = list(seq_dict.keys())
+
+ if len(name_util) == 0:
+ continue
+
+ # Box plot when Replications have been detected.
+ if any(int(name.split("rep_", 1)[1]) > 1 for name in name_util):
+ # Extract the values from dict
+ sorted_seq_opt = {}
+ # Number of replications
+ n_reps = meta_model.ExpDesign.n_replication
+
+ # Get the list of utility function names
+ # Handle if only one UtilityFunction is provided
+ if not isinstance(util_funcs, list):
+ util_funcs = [util_funcs]
+
+ for util in util_funcs:
+ sortedSeq = {}
+ # min number of runs available from reps
+ n_runs = min([seq_dict[f'{util}_rep_{i+1}'].shape[0]
+ for i in range(n_reps)])
+
+ for runIdx in range(n_runs):
+ values = []
+ for key in seq_dict.keys():
+ if util in key:
+ values.append(seq_dict[key][runIdx].mean())
+ sortedSeq['SeqItr_'+str(runIdx)] = np.array(values)
+ sorted_seq_opt[util] = sortedSeq
+
+ # BoxPlot
+ def draw_plot(data, labels, edge_color, fill_color, idx):
+ pos = labels - (4*idx-6)
+ bp = plt.boxplot(data, positions=pos, labels=labels,
+ patch_artist=True, sym='', widths=3)
+
+ ax.plot(pos, np.median(data, axis=0), lw=4, color=fill_color[idx])
+
+ elements = ['boxes', 'whiskers', 'fliers', 'means',
+ 'medians', 'caps']
+ for element in elements:
+ plt.setp(bp[element], color=edge_color[idx], alpha=0.6)
+
+ for patch in bp['boxes']:
+ patch.set(facecolor=fill_color[idx], alpha=0.6)
+
+ if meta_model.ExpDesign.n_new_samples != 1:
+ step1 = meta_model.ExpDesign.n_new_samples
+ step2 = 1
+ else:
+ step1 = 10
+ step2 = 10
+ edge_color = ['red', 'blue', 'green', 'black']
+ fill_color = ['tan', 'cyan', 'lightgreen', 'grey']
+ plot_label = plot
+ # Plot for different Utility Functions
+ for idx, util in enumerate(util_funcs):
+ all_errors = np.empty((n_reps, 0))
+
+ for key in list(sorted_seq_opt[util].keys()):
+ errors = sorted_seq_opt.get(util, {}).get(key)[:, None]
+ all_errors = np.hstack((all_errors, errors))
+
+ # Special cases for BME and KLD
+ if plot == 'KLD' or plot == 'BME':
+ # BME convergence if refBME is provided
+ if ref_BME_KLD is not None:
+ if plot == 'BME':
+ refValue = ref_BME_KLD[0]
+ plot_label = r'BME/BME$^{Ref.}$'
+ if plot == 'KLD':
+ refValue = ref_BME_KLD[1]
+ plot_label = '$D_{KL}[p(\\theta|y_*),p(\\theta)]'\
+ ' / D_{KL}^{Ref.}[p(\\theta|y_*), '\
+ 'p(\\theta)]$'
+
+ # Difference between BME/KLD and the ref. values
+ all_errors = np.divide(all_errors,
+ np.full((all_errors.shape),
+ refValue))
+
+ # Plot baseline for zero, i.e. no difference
+ plt.axhline(y=1.0, xmin=0, xmax=1, c='green',
+ ls='--', lw=2)
+
+ # Plot each UtilFuncs
+ labels = np.arange(n_init_samples, n_total_samples+1, step1)
+ draw_plot(all_errors[:, ::step2], labels, edge_color,
+ fill_color, idx)
+ # labels = np.array([10, 30, 50, 70, 90, 120, 150, 200])
+ # indices = [0, 20, 40, 60, 80, 110, 140, 190]
+ # draw_plot(all_errors[:, indices], labels, edge_color,
+ # fill_color, idx)
+
+ plt.xticks(labels)
+ # Set the major and minor locators
+ # ax.xaxis.set_major_locator(ticker.AutoLocator())
+ # ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
+ # ax.xaxis.grid(True, which='major', linestyle='-')
+ # ax.xaxis.grid(True, which='minor', linestyle='--')
+
+ # Shade
+ for center in labels[::2]:
+ ax.axvspan(center-8, center+8, alpha=0.1, color='grey')
+
+ # Legend
+ legend_elements = []
+ for idx, util in enumerate(util_funcs):
+ legend_elements.append(Patch(facecolor=fill_color[idx],
+ edgecolor=edge_color[idx],
+ label=util))
+ plt.legend(handles=legend_elements[::-1], loc='best')
+
+ if plot != 'BME' and plot != 'KLD':
+ plt.yscale('log')
+ plt.autoscale(True)
+ # ax.yaxis.set_minor_locator(ticker.LogLocator(numticks=999, subs="auto"))
+ ax.yaxis.grid(True, which='minor', linestyle='--')
+ plt.xlabel('\\# of training samples', fontsize=f_size)
+ plt.ylabel(plot_label, fontsize=f_size)
+ # plt.title(plot)
+ plt.xticks(fontsize=f_size)
+ plt.yticks(fontsize=f_size)
+
+ if save_fig:
+ # save the current figure
+ plot_name = plot.replace(' ', '_')
+ fig.savefig(
+ f'{newpath}/boxplot_solver_ishigami_{plot_name}.pdf',
+ bbox_inches='tight'
+ )
+ # Destroy the current plot
+ plt.clf()
+ return
+
+
+if __name__ == "__main__":
+ # Set variables
+ model_name = 'borehole'
+ solvers = ['BCS', 'FastARD', 'OMP', 'OLS']
+ path = f'/home/farid/bwSyncShare/Scientific_LH2/Promotion/dissertation/surrogate/data-borehole/'
+ f_size = 45
+
+ all_loo_errors = {}
+ all_valid_errors = {}
+ for solver in solvers:
+ # reading the data from the file
+ with open(f"{path}/{solver}/PCEModel_{model_name}.pkl", "rb") as input:
+ meta_model = joblib.load(input)
+
+ # Update name and Concatenate
+ all_valid_errors.update({key.replace('ALM', solver): value
+ for key, value in
+ meta_model.seqValidError.items()
+ })
+ all_loo_errors.update({key.replace('ALM', solver): value
+ for key, value in
+ meta_model.SeqModifiedLOO.items()
+ })
+ meta_model.seqValidError = all_valid_errors
+ meta_model.SeqModifiedLOO = all_loo_errors
+
+ # Plot box plot
+ plot_seq_design_diagnostics(meta_model, solvers)
diff --git a/examples/OHagan-function/data/valid_outputs.npy b/examples/OHagan-function/data/valid_outputs.npy
new file mode 100644
index 0000000000000000000000000000000000000000..6b0f28b72e05d71f6a0c7a677f416e9493301ac0
GIT binary patch
literal 928
zcmbWr{ZrEg0LO72hwflO!C*6B$mSk#5AqnsJ`@M&sf>sYhl<F9h%)A8J9bDdxG7w^
zy97i<41_CFh*9TC4ESnP2sn8dmJguLP;WF;&e1Gr*}u^HxA&`lOZH~jR(CgzTP3@&
zNL`>}zr|p$FOFvm80_M5)d5w$BClLkSoEK#=I>J%y>fL)zOv}mbHpN%fFb4z7*&k_
z8JXzV_U2%P7<+v#KfW1CL(siUbAN6yYPA#e($51RT~0_naFh$%(?50<HzXqdGo!hF
zyF<$5yrz+vV({*1x3ia?`s48))tSpIDk{Ir)qMO+0_%m+l_Sb9_~;z+$i@{oS}}Dl
zF)IWqB0sBm`}BbJ!)Hc&U;=!<AX#fpt;6~EHtPPlE{5{ONYecwr_@=^A9`LR0$r*1
zN+&bX{uiTfp5_ONE~#g-n}>!COR`A=85&FdxC;J6TwT`92%!k|=kk6US>nLNQ|o(<
zFo4Qmq<IJ-cE@d<Gqy!Q^=LA9p@IXeeZMgFfhSS}gjv36Tx2u*q6X?BVbctoNh84s
zsr&6#8ij`knz`z#&eg<8EqzvKj6~F-<k<^?K<LNKS5+=PBukC0S-o`RSPu5z$*>Hp
zZ|1MD?^%o5YZmX+A`y=NetSalvkwf?X1npA0IL`Qxf#M3EU5+uM-3q`B3d)n9ga_F
z^aNKs4>{yTzJo_W%*h=$n%)t>Ua0uwl8%R!k4)2^3JN&K(w0N5BBXsSORk>h;o^&a
zOPh|4GKYKrrqf&`pIIB~*A@w?T-9W2iN@GR@n3OcgV5@HQsOM<!i%|X?~KI*hZY42
zyMlzppw&;3s#!?C^~Ubw%vj9p&IYUyyAeB2EL6yC4EP@2Ok4PYPRP5weI7Qt!z1c>
z88@7X&b|ikX8&09T-|meEQwCMU(h*jc{>8r;)nC!1d4F4g@46R6@&BEf?rM3{y4J#
zPC<rR0?{$s@MKdc@j<xnw<%>D?3kMM>Dt7E!XzB+x*dltwcROqX>pjL-jHiMXvlr(
zypmVS0#QRU6EoI<o!-497gdP(zweGwOZ<p0#bU+OAQfao=aBB21SYSJ_9k8g_LZJ$
z*~Oy3fzIv;{Tig(UK$uNEO;naUDQW~LuO)Uj$UUXGw09DAHQdz!6i?0MGH~eo4q_d
O_c~g`67OrYp2WYvI)AVL
literal 0
HcmV?d00001
diff --git a/examples/OHagan-function/data/valid_samples.npy b/examples/OHagan-function/data/valid_samples.npy
new file mode 100644
index 0000000000000000000000000000000000000000..2a2a38e1a1718d94f39243fac3d0f047fdd870d1
GIT binary patch
literal 12128
zcmbVS_dnI|`?o4dOUr0@8`+Ymj9e#C5|XUaKv_{{B9cm_vI`|7J4CjK>)6WPdvk2g
zIrjdX@4xVU{r0+lxgW3BecjLNIj)ZyI_jD?=;<8kTm<f$*x$1gkmVMTH$5vL$}M1O
zZRcQT^ysd&-F=h)Tfb;zWp6@Tw>LMkF`>N+o|KRf<vuAbB+C7a`~UZ_#N`_E|9vuy
zvEr9acOMx?eW|r_>myWL!<u96qpPS+NEY01B_l8E3$afU@1UtA=ESoP5pdAbQ?D|w
z60A6{noqm-p|NBT-BK(WK0KFtVj(z$rt|9ZzB(%qe2Jmt$eBVc?0r36c3~16+3j!L
zNvwy1{ZI08y!x=|nQo6n*9tti9z4r9IS%{dqnbiQYH`y;S<<3o3cN#lejeoQfCG9T
zoW+P`z<hx31bNFGo;hH7O)p~z1&v-dmS1m0zh&*|BSVGAp%WPv6!#C>3;#vBi)_Ns
zL;Boyjzwr#&n^$Mn8Jj#Etx61bFtA;YSHfC3>cb_&v4J=pz;j^QL%%?=rMfK)ih`w
zADAb;kh(nt%Nv9hp6eahnXG25awrE5rF|0m!rh7|E1fl5^C)1t<8$Mg69bsddCRBc
zWeC{p7~3skJc%S`_B`v(IpiIbmV4`w3}JrMU0w}T?9&ePm~8qBmAZ;{qkPS1A<2JU
z#(N%m+G3QI6IStgUEnO6LM@co9CogO&$uNTEVrnTVQXyY(nYf-@P9`f{MWjTDA4Hp
zJLPyOum)-l9Fxq)To^GATkgg)TN^}FojzfAfoS6CD~<SB>?%9U1Ou_Py->e7_8+A5
z%-HX^F^S93hrURuGy{vshX(|{3^*Ts!NFm7Aq4%c9Ekosiwzd%LXOn-L-$_SlP>I|
zxb!omjAwirdYS%4`)81#j=`wXXj?T(beYdQyxIut<udPc=q7Qad(dHq_TQo``=o+|
z+Q7u0U^B=$3$L0-HXUcfF!K2Mm+3<t(7EW4_48&Ae6l{cHuR(!-tfP!+O5@#YucPW
z^5g|zlGFT0uSWu+M&`AszD01BasR-_H4VY0y(uz3NYFanzgJ9o0t$pJ8C7CtAdLA_
z&fDA>c(vi2c1dmp+=TBPONd$ryF}lqmIM9Rd?ayhID?GmQf|(y#jJtSvyY@fvkf@y
zuuY^fs0wW!>T#U+Yy@?~FU9Uho6$Y9j5<c%K;^GVF%oxraQOv)tIw`!w0IC(8dNn7
zP17U)9dR7Rl`LC3j^0Mh@>o$1q}0LhJuyGu`%uty^vUwq*F&JiaWYxvmOG4G&q&s7
z?}jeokH)X+R1i`y48Gz}4L)_xy51$sp;W!tb)C*3xYp3|TiK-_?waZji&9=;%eAZW
zjPykKnRe($-Sr`GDvk4}7FGfO#Ny#c;)C#Scu*x*Z4yJ{wQo7yUjnwIQSEzD(@^m7
zy3kurDzXW)&l&7n0x`7(vAMEN6mI_VI%O&c^PDeij?ztIR7UzU;h1`i4%*3O)>Q@9
z_lvwZ-a`R*r_bWHymZ9=U-G_=pZoB4?~i}8`Eht*#!2+5-Uy@<U)5R)Wgyd;+W(9i
z=P+1_Y^c#SjgJ~_IDKyxV{ZSwSY}NMJ~@;AD9U9VwHCr359}MpjN^Z7*GOx)@@KZ7
zR;CyDnLmI1sNDjY7uTX+wYK53K;4#G)nr(;^36yIZ-E)9k28L8b*S~c=a5(a0*0xN
zvq;nKpWn!fDB0?HXh@6Gzaz7iFnmspTrX3POrBPU#BSF@eU0e><%Cffs5Fj3;W3E(
zZ%{8!q!phIQw>;c8{lS()qX#@G4Och{&%Y`9g)qPtw(@)06uxtbsXs!hNFoJtny4$
z&`fIWaq(G$gC8s1^=?ie*9_kmV$Bf#;XSpPqCAcFZLFEo2xG`)Y<F2Dei~FpQ&Ogb
zhT-<*N<x|<87R+rbMVs|u$td%=iN%e+BWCX4uLVeuT$wy{IH6%e7DGQ_6@)|EaBeg
zHjaChz6JPG2T?d+b$xR(5DL=n|KVn&g1{c>kJU|WkiaWz^1h)RB%EFjxUP*r^qw@E
zqi?7{ccl5j*1>YPbYW+cY$6FS<dTJ~H|mk^b#PzNO@AoMw7#sM)d05E`}#J#7g3Q>
z*%dOo@wNHIMAlag==&>g$CKGaq@%oD(JvdpA=5(%+HU#4{rojw&twhKeLPosfZ74s
zTj%pby6K6+2^r7BTn5lP-&CR~vIv=H>Lbmgsi+nFZ9R)+4yq>~l^-9W;9%nyTX$PB
z&<9hvj)!*Rjg#kO!=hVY{|6<__k*Lbl|ka)W1ktc%`<HhxkpBd@Hl&l(+D(e<Bue9
zZ$P>BSF1xO#~}3S+OJ;gJP3IF+E~ba4ms=pbrWfQimp&NCbFRz=&o*ywYyk^be%7Z
z4_S`jeX-S_%%KaYtbDG2d(}Eh<nH^fn>d3v%LT=k@|JMu^4p>p59WZkW7ko$on)B(
zm7D%&HXc1xouBfH*1*xw>T@e4%P8;ZeLq%u5h5)AWW5=Q!y~LFrTdSP!D>)FJ>9VZ
zzqdVKdh5^yVNlMxFguCEKUnqQQ2{vavS&1A9K(+a<MgmF0{1plrC#h|AwCTM^0>Ez
zp1}4f1zDVT5DNmCU+5doLY^k;SZXL4crS{Hm^2RJ`0Ux_J$&>8VPeo$rjk{dmtqY*
zb+rpci7|~fEQ`pwp}3X9s0K&}x=Mn$XMtfPVpDu>6Qp<qo_>z6$Dh3NDR#MC;O*y_
zdZS<xKL76go*6uk%Bqh_Q(sWvml;DtcgH+%`s8U-n04@K$=v|5F5-!rGs6z9<KVo^
z@+)F!44JpElFcOtAY(eJgr92?@5WY4E2O96Zu|Y>d-7{B-Udu8cddX$MEQjqJ8R%;
zKC1`|PY?8ajyo?q=ELzeH^0~W$Dz^e-k}p7B-|NIFwC-@L%QBAIwCbixFr5)@|i^=
zc9p^(y1%Vh@V>>MSgQ`Q4o*yTwohW;?idG&kK?GR|C116wt(NqHV1@`m4er}T&A!(
z1*_DSU$OE`pw6bj^<(+VXxm`Zo53^ylMfpFWIQ_nZa=?pxS<)E(vN#&vXYR!!R1C`
zG96*Qe0k4;Umt{vy9!lLkbo`WIDzj~KC;C!kzC$)Lsj%k$`%$^==$&|RW)l6nj^p7
zE)FWjPU4}6U%P@}BXaB?$2)q0{t5Q|<)`~V_rDn31D<uTPj7xR*RvAunz-$6s2IeO
z&+Gdn9*~jmsa^mZ13fXlJW4YC{3IBL1svS+(iIz)RrXMm$jGr)bz(WMA3OUd-&Af)
z;bTS)qe8o3(DY7Y>SQ9L=wQtA$i}~*kVJRmnln8yUvE#U#!w9^9Xg=kV@f+`lMl=V
z8tovQ!vA=>@i#F0rj^~&9mkia<eQepi&1*p+%aFn4rGl9RDHOu7zQlE_5|G+#1)bF
zqi$&o1a8(d87`$0uqRT#&NiwWXE`GO(+({I$<K~OU&C`Tv-;UD!%NFZzU+~jlfQ<l
zMpI81SpK4s?#1(Lbj^UPh4rWJMdSXI+ar=@9e6Y^_|x}HMk4d}?Umw(^FZx{=xpqL
z3aD8;@Djro)bHtz;_6$5FUNoF|Btf)ww-GBUecI?@M9I#o>EkJ@+|4pzH7aB*y<`-
zOJ@SyeLF;>CPra@6idtQZF3;_hL9Cz77f9%`fR)czmd|MBu%g!#lsOT=Rdt^f&*ea
z_qTeJ@%TH#j;66$;Jkc((nh!r6`h=#hL>j{_^EE9_Q57lR(uylS?@rybYY&}zC@H2
zxx+X8Y8}fY;;z8Oc__~gZqE*0z-!eKbvrgP@V^2P?`FASU}8LJz^_kFbk^)qWBa^-
zuAg2#>W(0zpssdvn;j8vt-RJ2eb|V@?XnJf&iRm5l9{FAFo%q+yVRfmoWQCmhq0@A
zRLuR$?_pIq2Roy9>z9&7q5AmMB6sr{Y!}`Vtt~)LSiAN#bj4&5HEgfhWnLnqYvCo=
z1fLZ=caqoWwM{>sy^=AUMXL`!DMmlb14|(F!${?yZvhTnN%?!{{v7Q2Tw+|tM}gSw
zyk{Lh&jFd;04~`tWAs6Vi0Yg%bg*BRywBc=<;NMG3v0DOuCBMCNY4O<_>6RqO;Kt6
zvB*ugq8;l`eG^R6ZwF0&ORLZ`laRQ1r<E(P4@oE1uC4b_kc}tI_)^Xo6e)V8npV6)
zL9cpt*)5a!V^5e`5-%em?oiE{Qtn<%5OvQn8K1?=&Yz_>h3as>ztUy?>;F(a`9Uhn
z?RkiqKNhf=JO$D(c6SHZOhBjp+mE?F2hc4?DRpsY5q69Ed{qst!;MG&(wir!;H`CZ
zJFdmyy!vRoi`+cc73}TQ8z;fZ`+W`ENB%*TS;o-6BUB&~Wn?)t=isGIQ0rm<12Kdz
z{aBI4GKe0?_?Jgq#65|7?e{jeL9N70t<h99?tONA(oU%nt_{(z?K(q8lza19mhbKq
zX0vUOk3Sg%!B$h(7^P0scA|ff&#{cs|M_T*ru+ep1NMs{)0uGgWqEZWBMJMf_A^L1
ztV3q@vWd3jA~5X8+9}sP1JQ;b{Y|V^F+Y#-=sLL**?ZzN`6gy@IykPa*nA!*nr;P|
z5y(Kdec#Od+7LwR?e^Vh>B6fjG2%}KQ$f)t*V?OL6Q!@_>mE#B!O);LUk{z9BL<Y|
zU**ZC!2SHtz{FxYLitPWkNq)4FqzFlC{`TDY?6nmWcd`3q*DhUoApE5e@W|Jt5c{P
zlThG#(F4Djov5t2G69B1&owCRTR`_R=G%dHTQPXb+)Lhf7SHG2Uu%1{0A0!+a}K$5
z#Md{tFS)MJ_`ZWGo5##=5N(_OUb~Zum7#CWtVhkl>By7WzhwHs_s3}i^S_0-DHP0^
z-<AajWLF6AIuB1hTs|V8*$M1@qgyzg>(R!kY_>!r7p)ZS-UoiFL^oSMLh#3aaJKyC
zwca@fLgDSQW&wZEse6vR?bI}koY`9}cXR=ygPTdae*`1TWZS!Kyd%hJp6mOQhKt|q
z<%3KXsOU~GyLm@34;#kmFG=Jq;+uwPJ@M8GNDkC}Lisuk9u<9>CcL9y?tE*jt6?jo
zB|qbw-PQ+%^dD_!%L;*u7(4cWb|1?ROqRTh9m9t|m=0g3jNpe$Szir8Mu7XoeuuS<
zY`js~rWI?Qg|4gR(kbqYgqUQDkUNz%IDPoR@Z(bhC|Pz8hkw`Mxr*Yi3NeGAHl(}9
zG%OEHcNJ*8G#tR>&7DK94AwD@Rq&J88VMTpUDlQfb-0kOcvIc30_}=rDzaA>;i(5-
z>lN2_)bIMnZmZpaf^pN<W^1#^Z9a)gguk#{^5o&|z#&Y_4Y&qst;qY_VsFy(WhgEk
zl*GnmO!STT_NcxXj#NGm`e`_UwLAEatDRZ`g%m@+oaA{}VE+@&5SEBlzYSD_q9`cE
z>MuSNSdVS$bAIuiYd9$U?b<F5CW1sx_sfEcDXdml<7B%@>#J6mZr){G05Ks(=2Ync
zVB3}A-Cy4WTQb=0t~ve$hx07di)CF9uAORB>pqLUmZolk!Q<GfG1$IMClmcY%&r#v
zw}4-&tlS*DLoq!h@NmVsDumfHE#VSHVEA8x@6)f{=(Rm*5@@_PclX)h!W)@*RCeXn
zu0KtvSaDU!z<(7kU8?os{xbw6S*0C20>)r}AE~?g>I|k8a=pk>jR)3loo3szWk};1
z?QxYF1|IcKif;!g;H<W7Q|0G6UVYc$S#Ma43tnooCT(*_s<WN8;iF>VnDuQ}&3QDp
z{Fi*3lb+~7@IM_cMaDV(8>uB1D!{d9uj{*kMo{I}BMD|V!t}Ar^W|~FVBfR*)SlU4
za5p&~J9TCb<M&$({GHkWi?KDj__|iq<_PIZGc16cN5!A04wT?jO{=tNYzy*UGW}^<
zL*sc$o~{d9L-6eJ8)r{IJ65tjsYw-@hf$&HOkPj_BAvk51iz_%w7)y(p83Qa&S%*F
zOZ!TK?6!ti{GKCF78o~mM}7wOve<OIiz>%T1M`r^17)bm{AO!MN-IdX9HX>4Zy_Xm
z7t>Sy7t!oGANBILH7HU)(UW9Z1YCUwi1UH-xKq&nwvXT_b{eUsh>nvmVz-rBzSAru
zo%qi0Dct~aTZMM*E@{C4Y2GLEoa@N`qR@F#b{>Y-t+aHjMo{8ri|EgtlOP(jz*y7O
zkEK>ef6hzyBfUjkK<c3(_^?#+MS1TCU}z%$QuZ)B)o$9~XHkr+1H1XAN@h_1r;AjI
z>>$#KyvR51=>rb2Y_p}E?eN_5O0$<C6M^mQ`-K+qDp)3$tvBCU#;p5k?Hl{1vHqw0
z$#bjUurH8n?vo@5iG@sS!=J-&58rJCBj#=(tT-?Q=Zs+Sv$?+^W780I>Y$0YNI$l5
zG#41XWgudD(vadv57ZqC)~}{kf;O8_<k1}iaOc27I>sY&SZ+!R-M41~lxw`I=ns&v
zeU!ZU{`e%)2fw+H8MOe`gU3?`>qn4#*ZW>C#xU4ba!&r%g%u2<`*8IxOFe{qGn$hA
zIFIjiWGNCdHQ?=RYNGUV6GsQjO4g{u*qOldK(A;X&UnU~vUH8Xs@j`?5R-mT)aDU$
zcsGxhDKCVNV;SVBCVC&cz(7=I|M#|4bQpXadnJuOG{c>{T-!LO$$02F|I_$yzaY80
zHfe9q5{BMy`oT%x2c~<a*5}IyL3*#F8IvIyIt3OOwk_Ag@a@`zlIOM(o^ubel}UDk
z>}u6F4qcj_Xx93pE4hS5I{NiRoReT%cb`YqX$G6F9DS61s|C-rn2nvcodk2^V0{N>
z68Kxb9(^3#kN+Naa+VdQp+$yXfGsH*sykgCj(sPBwCwWk&i6E4RV(6r?CTt))lDij
zeqG1lfq_qsHy7dis*fB0>nXU(hBKCOtC*tr#+F590}lVdZ+d)vFvR5buHc;yK7N0|
z!A`jwSx6N!X6Cb??|WE6(7zKrnu=KW?5V;%l^?UWpIQJb?a`&JVFO4NmmM4MprGyZ
z+I87pI$}7w{`jo6gePvBTigwrNBNM{eOd#vpt~<DZ8BsLLwGLd7V=F1gD@kz-9R7e
z{cu0^;AtoDNH1G#?^%O|CpM-~*$%I&t_%KZCSu}ClJ`f31vqEPb*>s|{>QVpoWif`
zP#EH?@v+SdnI|UiSU&7W5_t<Bm2m)Cy!X2IByWO}Qp)y|pW6_(x6ro<41;HfpEN;d
z2-i4_{~FyM1(!kPZkL^Xcr4qzDb{`*TsFvVB+fP<NM&)E)Ak#u=<{1s>rv~QfNao}
zIUq%SX6U=v4{oKQM|WEkg55d!Esd$`AhKmmEt9<*AKq248t|FIpwWwKU#SyV%2afx
zqvkKf_|2$2C>w(9&n9)BYR+PSoKNhn{6%nC*_k3-NKd4f>?8d#o`f&44l>!c1K4ZM
z@OAH%HjrlZ3A-4W3!SfX&5S+fL9;NIF}KSR0ym3O)qd5((c)n~d+IWJ#Z>Ss7dGI@
zg$Bg~(1Z8-*WRv_FF^VEyd3h?C8X!&JROk9Ol&UwE#`Q89)(2TdoOsk;1{_m-E*}x
zJ(}VPtjyz>Y`RYPccK_&T7S72Jg&ydV5J=yJyc*UJw3?xmX6SIQSH&AmPU-e)%jB5
zeIKX@Bz24b8%2%qh1ZX=S5f&6>3VZbHAJ}DmYdFWgG`;;f-x@%9@Ibo>an<nHpI{h
zs~|EQWF6O1R4;^<fjU*y+(B&N=hM}8AH=tR*IHbK*I>H`tJtq!vsmstS#jgZ64X5s
zv=X<Zf*h~iMU!ihpeK~<aw=g4ITUY&d7NIvS=Z&SY1Dq0au#5FO`5~v&M7*vfJVq^
zr$<4)2As;G*D1Kw398*%rws1)Vbrxu>#ItOFq0fDlY7Y#bB3)gog)UIVk)9rI(-_N
z)ax7vSm*Igqkpj0CKGY?&qc==xk~I67x0gpOhAF~oZCYhBs}veET?vc3gcov!sN$k
zSkK*LQtU_bY0~}#P!)RMp!u@#JL?g+y#2DS2hG>F60DNHDw7Srghu)cG<$#=RjEB3
zOhs0<Bbo$-c4!cMSL@uejDFw7e;=U^z-q{qlg^g(gdMSJr5TJ(XvF&E8n0gjrZPQ`
za|xP;j)#X-iY!WCP`yFm(0T$Iwy8v=v`(R4@>z$`^m%mPy!bUFyaP5r<=(woH-R&y
zH&NvI1RS!-%aUbigd2kC*Buv!ATH+2$+I@A;P;zIac7vtCk+y_PFdCX?TNCV<wiO3
z9Fe|Y8`g-M{JK)LlEWBsXr|%-?Yn>H1#+52_M(yTjT;wqGqA~R^l9#|Zm|2QcY3F{
zG`yQJc~!MSMW@MNhTW@45Hefz`SGg^_z}OOO;olHs9#Mz%2sQEeXVbzhkpYpm1p@&
zXGoaQ${S^{x`^YEKU7YzCc&~~|3faiMo>AjWB%CiCO)@U?J(0FLB>vB<67xPaP6o`
z+#554oPS^M{BN}%2m0dfobjB54~-r!mnA2`QqZbw=d)q#(ATU8_>TlN(lR2af0ZKF
zMt|~jc{7xlhaCwBT7|s$+*PN7UdS#=t$gV*f=YfWf#fT7Xl|RCx8?INWao^e#|;p%
zVnj5HiGKmE<cAh?=eI-Yb}dUDohBR&9x(hv^Fu{B-*oty(-Y2Z@UG9%{9A1axtk&X
z4dJg{cP{w2_JG^<dG@TZ5gOi>wscYQ@i1jIY0`fP=c~^aRZ_c9*PJv?bfx7pT8`PK
zU;0ob=P)0=d@IHqlrpt{Cu4qMH-Vn!TSvDZixg}o!4QLOPHNczdaZbAW?R*v@wPS>
ziq|YSxc~QtS~7>nc5zQutHr|W%eG(KH%nlPs>FehJBN{;)wouDKM6A8qSFLKXK^@@
zt>=A71@51+?6o5fA%o$7X840K42@2>YL=1-a=xsIXWgl=8QEV=V3`5xr_C-k_i=DF
zVq7kK*aSzp4G*7<?8lq>vbA|vC@A0+ds>jK8krjX#kjXppboo>eal85RBzXX_2Z*7
zJk!1t7%%~BPxf&5eWW9t-TsGA=r#zfLW`lThuYxv?8}MeXRYvn{r2{kPRrPA{Wqhj
zy8^O3|9shgp#kNM+ZAF{(qY?q=^p}2Z(#qj+l^+827GdThs$;vp1fP~p~`VBzzxg%
z`QAhA_|iZ<;($#*ysm$(pZ8!M7S|7ZY50=ikfHVADho0w#TN-nvAaR0oozz7z#@kK
zmX@;QECk)}b-ttKt5C2%ZmNv87aL~#wwDf4p{e_Qrtj|(l+keQQ+Uyb%=+XJ!#Wa(
zzW;ui-g_9eQtD&OLdW5nqKu5tjbQ8(C}KLZkcy>m!w>t;zk|m|8Pgh#yD`$FA%f>_
z39e?JfBgDlDI}k{TK(VC8Jcd55YcI_2Ldrv;?da_AnkWN4un}eZ9_GBRY%3JT+{4;
z+npf1(<rjaw;AiJ`?jZ1M)8ltr~oSV<0sQ>Ndfj%WGh@Ih?A<&>QwW;<hB;vCN_WE
zfpZnowrC{?3{$}Sc*$A)`X;co7&5uQ+=hXZGW&i<Qc!;7Z?En3D(rlS+Wjem5Th^N
ze|Nqd+~1uPOnOMm8Q<zuse1nfA1?=T&rmCnhW{P56YE2hjwIjav0CJ_ImQuA<Bx>_
z@5Z{D$1puJb&6+b4ft;JxPG`(3gL@4)m>iHfeAO)Snt+8c)=s3ct5lSyR{3VHmArK
zX;<winbL@y+{O)Z8Kd~@@|Mb|zDBSq>TkEDOhfZ--FBw+G34%eqV_O)10x4T({wkp
zA?krCmvrI?mU9LN^7L(D8r#KYDf(9Ua_N_7N5&vFT!xOGv_;GcnOYZ*8inG`4EyN)
zKcPjc+<q!?3b!`=@-jS4!OWKR51`i$;qp3p*guN=Qh)s;?hm8gXYH&yS{(>*o^*Ta
zSB4utO$RRyk#W!L?S!-2hwx5QMDU9-3hXfqw;7}~!WnyJ`6n!WxKphFzqrppcZmj<
z2rai{uGva=#c31%7T*KPwR(_`bh34N+=uPQ2&_W#WRzvII)S^)At^xkdzaH1`Yrw^
z%i~Y$yV(c3BgE&R`PlJGx34k~9&v(mN@E{}`(_wy<)L6MY0YteNEM#1qFi#YZ3b@M
zEh37n9Y|N%aL>-a4>=56Q28$%@yfp8@B-U*fJjZf@a8G(yTbLDv~Ls2rD{}AlML5$
z)ERP5FGI?%cNG~1W3aYnYaG(ngRy58hWAh?*zNX;ZJ~swJLK6|33pq-ukq5+v&OXi
zt(ZtSc&7q(N=B8;56^+G>A@(KZDcHB_+}t<X&(58)D1qZPt)?6*LHNAlW1xm+q6Hm
z65gG+&{oeHN2}?$9aaVtILO-bA32$f1(y>asTC0+;cU^=OX+r0-cpe_c$bbyag$Cx
zxiO6UW9v@wI#5vIAA^<5l_8jYGjQhK!(l9GADvN3BjdYA_b(XCY+#eU#_SA@H<(=0
z3|#W=1|IRev(ptMOlr$XZ`B}Uc5LHkQ-cAt3H>gy?6Cn!NblmKGzA5Z?daYQ4<OqS
zsg**?47h1M%5WY3g$wM}9YoIU1V@nx`D@o_q0!OolbygixJM<Myfm%HJ&aXjj>2Sg
zx*YXb&zg>yMO05bB1G#KDpm*EB%2`M;TD@#-xM0pcZ+IiOoER-*L7kX)A7XQmECgB
zJMo%K(IH=2PWg3Zq&2W%6w;3Fdo{eyKqOVn?vPe(1lE0Pk#oWAa8Pp2z&~jNS$>D7
zX=u=Nl6Tk}we>}mNWb(u&1wk_j8SE|Uy+d8O;zHe&m?*uXJ^SDj-t&e=>29KXv0?l
zD&*V0{y>JQstbc~A_jivzq77ei3LKL{)a`W7?E)`HO82Z5ayep{DYw!*R)!+=ouT~
zUbax$i1-NFJDbQUUl;-t;j*{0FNX1o$mQ63ejUg}Yzu4)EJd$y|JSK%<)~!ocxTv;
ziVcUBf2qd|Baz6V_TyzOC^b~hym2GJ^z*ENh1*N;?AkV}RAL*lybXE&EutHqS~ux9
zi?qYpB*&Mm>SJ)+?t;eI+yvnHd3AD-o&+sBzH#rSn}BCJ8+Rk#N21`2^fz6FGC1V-
zC8mhjfPo|dy}yUsAo@nuJIVYiU@ta25^$LW^~^oLFVrjm!(e&Dko_Dy)??cnPi=vr
z@GHCRcNM{dy}U2xKIOnk<`V&mxP+L#zf;$C3_TBgFg<*I16&tgE5_I@qQ9$I>0ho5
zyy#|QYP&d(TAuSv2JM5e_C(!MQ*#p4!W@n^U2KE5$|(mjg;&sQ>d1#3K3ynuV%6N{
z(LAO)DAfGeGXUDmaTD>0CEyePe&NCUENnX7Z`Z_2#g2>RhJs0hD6RcerZ}+<yl`z-
zWgIP+{m7BtFFOeP%zkfAB#pxYSDDJKJ%4CE3s0ls$P%1AoO9R2ZyulBm5KK1ZNY|h
zAxE~ZF_a!z$QPm?hBvQFpTD|8fl(6{t4d}H7{%MvuP3(Si4AFo^3hHtuB)lO%3~th
zHU2GX-N8U0>wVAKV&0BM4>P~EP<nycicMoOW)6(cM6mQ%*29gW;Iq%mM^Mc6$d@A5
ze&n=O`LVwy8;r&sxIXThhd!>Xj&IeBgzc{JFWPtYz__PV^rfye(5Y>;Yf-6&?Ed%n
zG!9omU-{wKWWhmHkbA?u^t~VN`xO_3j4#2-z<|QPjmwbyR<<`phQ^Ol=bj`OwV^hN
z+h=2H0XizuCX{vR@u5xw$l4NNsxVW{-(?zp>KqBY<UBx|bE(y^{<Hv=Qr?H%r#tZ5
zHC;b5<`HyCzjlaR)P&WI0?)6#Sb@L=Wo4l8r*Cf$^~WC@Lek~1n@3$o;j7ZN$8PV*
zxL_am$tP_U?wyV+{36kVzPD{IFnwvm&+(;hmu^t7z<2DC7*`ug-4=MsB3Fi?j30~5
z3k!kegv1qv=zehFVsv#6>qm}{N!x0b7zmV$Bk#T@je{>qPW<J6T^L2p{oyITjd($(
zT;}aMBjK6NzAbBX14xKEI{Rf~3~V)3c?Dihz?~tZJ>Or9U`VF5px)8~4tdW0xpKY@
zv{SQeWEX3oCePH@sbT=dUW}YI_F4h)8!=b5=Z|1z^1h=p@%>;vA!{yko{XJ-LH`2Y
zmmvKG?eJ^af3fVwXT?V{8R$#8!%}{xAC=shPkkQhfdOW%T9bo}L^iJGeeD%v_`_&<
zyKG7yM!fLl(4L{7mGqU}H<#w%(c(TO6E9jlY6*6yy4PTToW142^8=tZAM!{_bOKjh
zlImZC&tRF##|v3!20-y!$KxG4dw~)f96j9GN7GHUmXFN(ATjp4++JGFrz5xj^zM%%
zxW&)gMJ;C*lNT=ulo0xX`3E0|_x4E|pE@{H5HW&#<*gO6Pfy_!4!Zj%FZJO`iJ8bK
zhH41fx_@_E?Iau)Rez}%Itx?!`rp2c51>Qye13CyGjO_p(hU}C1MimhyF=|`cxz^<
zd$V@{xayWf6l3f0&H|C`ZayQS?wOOs(UYs#r2kvRBfkU$JXTM}UzviQid#2+)_H=V
z_@D>Yu0e8I+E20<8GX70p8xyTfc~FvImj}NVt~x@_oMoA&~8G`&Qu;pi@Tjq#U{wW
zl)Wav+d7D9HJ5B=BRhcet%aS^{4C_yY<JIXqrk<(o|<~5sqp#XR?dxsRCwVj<+je=
z3)9@1SLKRmco3dPnA=!_DAjekUxZ#fSn!=in}#o=cQoD(DlcQk)yJ(Hw^yJ(Vg2+I
z#wNILn5x$&iQsrK<3*PU?Oa@%#4ac<flZUB)U$#mnD9KHK+Nld)33P}^&gfZw{ZvY
zgJL!u<vI519L?|kl<{ungvJp3)E}Gmv=~G_t@^qu@n)>ub=&ls+XR@T*{)P*RzVbZ
zr@@CbCSq?^c8l`fL4d!PcTkU!!KA#cK$t%sLf+VW7Bw%x#=qc+|Gsr#(RH@Hzi&)q
z_MZeg9x#S)6J+s|r)F@<h92Z!H{$|>UO}Gz7%J{H+)$_Ge9;eheRN!>z~w&g$4h#3
z$RBn5-618Kk5$kbG;$O1fpQL?qChI%d@XwEcf%$;myD_oELcVh2wQZz+>DYD@~?zz
zhH<0__&*Gi!DW{#)5&foLO8>RK!M6GXj6S~SZHAdo(5uy?dd`A_IJNHMXNLJ%55t-
zElfo3*S1Rcs)w-Rpjx4;RzGxpGgsvIN=34k``(zN<1q5(%}$FCOvG*zGxghwB)Gmq
zDY?I66m~Pe^f+Kw1`mGZy*<Ix44ty-m*XEzL)&rrg!NquI4>b5R>R(l1m;G5<AxSk
z(KNkp(peAM`DIoTduMQ8*9UF7&`J!|-!dAkItja0oOh(p&7olCq=2S99l`$Qz?$jB
zDyY|(SKCG7`}#rwr`<UQvFUVtMOk1Y-fR0(DaKs_-3Q-g<i&M^^-oWV>Si-W7U=nr
zOQ_JH6^w7>`@kS-K2C?02Y))fJ379Y3}=yPm7KMJPgxqWMG9)sPyB;n#^?eZ6!NzC
z`DX=fBN{lgb5^jSLq%dFoC<26GOSE@bs;Z$i-mro02lM0fh|9}@Y6(Z{j-znXxSlj
zQfG1$Pn^B2u+O^%`c>C+UdMJ~`cut_C!F(eH8{@cUhf3%7;b!`V>bv9^|=iiiappW
zCb>09YZentSH!Cf^6>oA42R1K)4=1YuT^<~rn|lj84X=e#&cBLQcnFr813g#U!mQb
zAPp7yD4r>-doIIH;h0627RR>f8-q~zI6R>6d=?UQ*qVh+7eGa+ox1U!3I*!YeH?oe
z;kzp5ds0Xv3b!3m*8aT=)*M&F#f2ttBizxY-C_aR*ZVzNttU`uY)_3>T?@#tQ9A4l
zn$gfz;&hJIBCKtSd8WoM!lda%l~tM#q<gt+pZiKIBspw;qn^qJnd*uhtG+pWle}A<
z;d&QrKYxc;d>;u5A8MsP)%l0Kvll-mdGtWl?tSx{Z~NiTQq<U_#R?4UkNX(BZ65XC
zJ&sqSj|8HYX$fa@Ip!UFF88Kp044}}3-OVo;Bn!`6z|X^`kgrcC{CJzNcZC9_CI-4
zV9CFB$+pA?qub&o>|V~JDd&K%|BpFv^kgzp(3r(vR|3*w<x*h#E3>ABF%k;gwKic*
zh=C32zcW#cgOHkk<6d<25_WpN6k`~&M~zw2Q1)>u4v2ho_>kL;5f8iyw;L1jpK$vy
z(;^YRf4FFKo*E5XPtOx0yGrl}$B{Gl{zbw0M3cAIRl?y??cka43y5#H#`iKgEg^U3
z-9Bc`K4>zf$cfrcK#A+=x|=Vj;19$5y-5yzDE{M5aDF@$+5{Gzl=@cDz$xtfHyYk*
ziFES?N^SrnW8H712i-XRd-yEFbO%&^G)!|m-vbM`ns}-u2e7exH&NV&iBM!v4N`l)
z!!uRRS|(EpTptqAJ+?f6ds$zeasOHk9D5WT?2<cyUVXEd$Ttkd39Y&LA!J-Te%J46
zeIz!0{^@ySj0)dx(LZz7pkm3M(^?saDKO{QY@pZLfZq-$5KfBCL9zd-EsV4|2!YTA
z_a**L5PK73+fqIXhs@3Dbf#<3CGhJFi{M#oifg!=^mqVvchAdxaE--tKkSY!E|O6`
zma*$`;SjbHCl0fwj=(41Q(8~Y5i#wIr_zp{b--#f8lV+P#c(y*l%#8oQ0~dZxHUK(
zx_)up{>?LuD$dRN3Dh=pl`>adu<gg`^0CDWeM=xZw^OKf?j!D;6;Ii4tiayCc13J6
z<tWVcp6Ol@6LBR)^g@FC6#TsT^Y~EPFBBxrYVz`s!TGRUto6}3{MS^mGBMBzHWl%9
z4K>T4%x}o7*g!{oZ?fF6&3_W#7<RDSlcK;Y_Xf&WCkhHWSTvm)?gP6+%@^CcMloS(
z*Y-TNZoFn)Yn*;_7+mW7gL<1rAgJHm>HTm96un*Oe(u)^6%i+`jyQH=pu&&VTJ{Mj
z+AR3YV7Ly7ry^Dbe3md@t5w+NKr<W?P><+zsDt}(Am5Q;88qY5N8Zuq4ZfAz?vl<K
z#Df`M(tBmbU_gk>%n(2Uk=(IAV!<<*`^ifF##j?R5!ja+b*vSSr4|QRGg8oZM(VZS
zJdJlglM9__SOVenN8!mX<)Hsq=>)fY1vUtUJ+xNoMOi;R#T2s^v=+ADc`iNx(vQ!^
zjM98!ew)8ip;EN@p-h$|d%w=$vXO#=Wa$#xejppy`>#N(78~&_2W?)|yO^}fl?oCH
ze)NS=R5ZG)_?Cr(irEVB><U$*Ku*+m6a2M~%kF&^1!?POeK$8sZ)OxP-jycUBoAOr
zE`vs)$2=xzY=kr`PJ?UWez%|0Rh$%PyqYgm3S}LR)bIwx3vXwvoDx>BaC7?7z|=N^
j>3<@p6f7r@zw6=9k<NA4UA1LdL}Uo;f4(8c(NXY!$*oS(
literal 0
HcmV?d00001
diff --git a/examples/OHagan-function/test_OHagan.py b/examples/OHagan-function/test_OHagan.py
new file mode 100644
index 000000000..4715ecc52
--- /dev/null
+++ b/examples/OHagan-function/test_OHagan.py
@@ -0,0 +1,197 @@
+#!/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 = 'OHagan'
+ Model.name = 'OHagan'
+
+ Model.Output.names = ['Z']
+
+ # =====================================================
+ # ========= PROBABILISTIC INPUT MODEL ==============
+ # =====================================================
+ Inputs = Input()
+
+ n_dim = 15
+
+ for i in range(n_dim):
+ Inputs.add_marginals()
+ Inputs.Marginals[i].name = "$\\theta_{"+str(i+1)+"}$"
+ Inputs.Marginals[i].dist_type = 'normal'
+ Inputs.Marginals[i].parameters = [0, 1]
+
+ # =====================================================
+ # ====== 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 = 'FastARD'
+
+ # 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 = 7
+
+ # q-quasi-norm 0<q<1 (default=1)
+ MetaModelOpts.pce_q_norm = 0.5
+
+ # 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 = 100
+
+ # 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 = 'user'
+
+ # Provide the experimental design object with a hdf5 file
+ MetaModelOpts.ExpDesign.hdf5_file = 'ExpDesign_OHagan_orig.hdf5'
+
+ # Sequential experimental design (needed only for sequential ExpDesign)
+ MetaModelOpts.ExpDesign.n_new_samples = 1
+ MetaModelOpts.ExpDesign.n_max_samples = 500 # 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 = 10000
+ 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 = {
+ 'Z': 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)
+
+ # 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')
--
GitLab