diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index 820bceede74ff8d64957fc7dc125858f1e99ffae..a2ff59d7721b5e110122548e81c164ab7e6f0641 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -1,14 +1,14 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
-This test shows a surrogate-assisted Bayesian calibration of a time dependent
+This test shows a surrogate-assisted Bayesian calibration of a time dependent
analytical function.
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/
+Institute for Modelling Hydraulic and Environmental Systems (IWS), University
+of Stuttgart, www.iws.uni-stuttgart.de/lh2/
Pfaffenwaldring 61
70569 Stuttgart
@@ -18,149 +18,149 @@ Created on Fri Aug 9 2019
import numpy as np
import pandas as pd
-import sys, joblib
+import sys
+import joblib
import matplotlib
matplotlib.use('agg')
-sys.path.insert(0,'./../../../BayesValidRox')
+sys.path.insert(0, './../../../BayesValidRox')
try:
import cPickle as pickle
except ModuleNotFoundError:
import _pickle as pickle
-from PyLink.FuncForwardModel import FuncForwardModel
+# from PyLink.FuncForwardModel import FuncForwardModel
+from PyLink.PyLinkForwardModel_v2 import PyLinkForwardModel
from surrogate_models.Input import Input
from surrogate_models.surrogate_models import Metamodel
from PostProcessing.PostProcessing import PostProcessing
from BayesInference.BayesInference import BayesInference, Discrepancy
if __name__ == "__main__":
-
+
# Number of parameters
- ndim = 2 # 2, 10
-
- #=====================================================
- #============= COMPUTATIONAL MODEL ================
- #=====================================================
- Model = FuncForwardModel()
-
- Model.Type = 'Function'
- Model.pyFile = 'AnalyticalFunction'
- Model.Name = 'AnalyticFunc'
-
- Model.Output.Names = ['Z']
-
-
+ ndim = 2 # 2, 10
+
+ # =====================================================
+ # ============= COMPUTATIONAL MODEL ================
+ # =====================================================
+ Model = PyLinkForwardModel() # FuncForwardModel()
+
+ Model.link_type = 'Function'
+ Model.py_file = 'AnalyticalFunction'
+ Model.name = 'AnalyticFunc'
+
+ Model.Output.names = ['Z']
+
# For Bayesian inversion synthetic data with X=[0,0]
+ Model.Observations = {}
Model.Observations['Time [s]'] = np.arange(0, 10, 1.) / 9
- Model.Observations['Z'] = np.repeat([2.],10)
-
+ Model.Observations['Z'] = np.repeat([2.], 10)
+
# For Checking with the MonteCarlo refrence
+ Model.MCReference = {}
Model.MCReference['Time [s]'] = np.arange(0, 10, 1.) / 9
- Model.MCReference['mean'] = np.load("data/mean_"+str(ndim)+".npy")
- Model.MCReference['std'] = np.load("data/std_"+str(ndim)+".npy")
+ Model.MCReference['mean'] = np.load(f"data/mean_{ndim}.npy")
+ Model.MCReference['std'] = np.load(f"data/std_{ndim}.npy")
- #=====================================================
- #========= PROBABILISTIC INPUT MODEL ==============
- #=====================================================
- # Define the uncertain parameters with their mean and
+ # =====================================================
+ # ========= PROBABILISTIC INPUT MODEL ==============
+ # =====================================================
+ # Define the uncertain parameters with their mean and
# standard deviation
Inputs = Input()
-
+
# Assuming dependent input variables
# Inputs.Rosenblatt = True
-
+
# for i in range(ndim):
# Inputs.addMarginals()
# Inputs.Marginals[i].Name = '$X_{%s}$'%(i+1)
# Inputs.Marginals[i].DistType = 'unif'
# Inputs.Marginals[i].Parameters = [-5, 5]
-
+
# arbitrary polynomial chaos
inputParams = np.load('data/InputParameters_{}.npy'.format(ndim))
for i in range(ndim):
Inputs.addMarginals()
- Inputs.Marginals[i].Name = '$X_{%s}$'%(i+1)
- Inputs.Marginals[i].InputValues = inputParams[:,i]
-
- #=====================================================
- #========== DEFINITION OF THE METAMODEL ============
- #=====================================================
+ Inputs.Marginals[i].Name = f'$X_{i+1}$'
+ Inputs.Marginals[i].InputValues = inputParams[:, i]
+
+ # =====================================================
+ # ========== DEFINITION OF THE METAMODEL ============
+ # =====================================================
MetaModelOpts = Metamodel(Inputs)
-
+
# Select if you want to preserve the spatial/temporal depencencies
# MetaModelOpts.DimRedMethod = 'PCA'
# MetaModelOpts.varPCAThreshold = 100 #99.999
-
+
# Select your metamodel method
# 1) PCE (Polynomial Chaos Expansion) 2) GPE (Gaussian Process Emulator)
- # 3) PCEKriging (PCE + GPE)
MetaModelOpts.metaModel = 'PCE'
-
+
# ------------------------------------------------
# ------------- PCE Specification ----------------
# ------------------------------------------------
- # Select the sparse least-square minimization method for
+ # Select the sparse least-square minimization method for
# the PCE coefficients calculation:
- # 1)OLS: Ordinary Least Square 2)BRR: Bayesian Ridge Regression
+ # 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)VBL: Variational Bayesian Learning
- # 7)EBL: Emperical Bayesian Learning
+ # 7)EBL: Emperical Bayesian Learning
MetaModelOpts.RegMethod = 'FastARD'
# 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
+ # error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
- MetaModelOpts.MinPceDegree = 6#12
- MetaModelOpts.MaxPceDegree = 6#12 #12
+ MetaModelOpts.MinPceDegree = 6 # 12
+ MetaModelOpts.MaxPceDegree = 6 # 12
# q-quasi-norm 0<q<1 (default=1)
- MetaModelOpts.q = 0.75 if ndim<5 else 0.5
+ MetaModelOpts.q = 0.75 if ndim < 5 else 0.5
# MetaModelOpts.q = np.linspace(0.25,0.75, 3)
-
+
# Print summary of the regression results
# MetaModelOpts.DisplayFlag = True
-
+
# ------------------------------------------------
# ------ Experimental Design Configuration -------
# ------------------------------------------------
- # Generate an experimental design of size NrExpDesign based on a latin
- # hypercube sampling of the input model or user-defined values of X and/or Y:
MetaModelOpts.addExpDesign()
-
+
# One-shot (normal) or Sequential Adaptive (sequential) Design
- MetaModelOpts.ExpDesign.Method = 'normal'
- NrofInitSamples = 50 #50 5*ndim
+ MetaModelOpts.ExpDesign.Method = 'sequential'
+ NrofInitSamples = 5*ndim # 50 5*ndim
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
-
+
# Sampling methods
- # 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) chebyshev(FT)
- # 7) korobov 8) grid(FT) 9) nested_grid(FT) 10)user
+ # 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.SamplingMethod = 'latin_hypercube'
-
+
# Provide the experimental design object with a hdf5 file
# MetaModelOpts.ExpDesign.hdf5File = 'ExpDesign_AnalyticFunc.hdf5'
-
+
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
- MetaModelOpts.ExpDesign.MaxNSamples = 75 #150
+ MetaModelOpts.ExpDesign.MaxNSamples = 20 # 150
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-16
-
+
# Defining the measurement error, if it's known a priori
obsData = pd.DataFrame(Model.Observations, columns=Model.Output.Names)
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
DiscrepancyOpts.Parameters = obsData**2
MetaModelOpts.Discrepancy = DiscrepancyOpts
-
+
# Plot the posterior snapshots for SeqDesign
MetaModelOpts.ExpDesign.PostSnapshot = False
MetaModelOpts.ExpDesign.stepSnapshot = 1
MetaModelOpts.ExpDesign.MAP = [0] * ndim
- MetaModelOpts.ExpDesign.parNames = ['$X_{%s}$'%(i+1) for i in range(ndim)]
-
+ MetaModelOpts.ExpDesign.parNames = [f'$X_{i+1}$' for i in range(ndim)]
+
# ------------------------------------------------
# ------- Sequential Design configuration --------
# ------------------------------------------------
@@ -169,112 +169,115 @@ if __name__ == "__main__":
MetaModelOpts.ExpDesign.TradeOffScheme = 'None'
# MetaModelOpts.ExpDesign.nReprications = 20
# -------- Exploration ------
- #1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'LOOCV' 5)'dual annealing'
+ # 1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'LOOCV' 5)'dual annealing'
MetaModelOpts.ExpDesign.ExploreMethod = 'latin_hypercube'
-
+
# Use when 'dual annealing' chosen
MetaModelOpts.ExpDesign.MaxFunItr = 200
-
+
# Use when 'Voronoi' or 'random' or 'latin_hypercube' chosen
- MetaModelOpts.ExpDesign.NCandidate = 5000
- MetaModelOpts.ExpDesign.NrofCandGroups = 4 #8
-
+ MetaModelOpts.ExpDesign.NCandidate = 5000
+ MetaModelOpts.ExpDesign.NrofCandGroups = 4 # 8
+
# -------- Exploitation ------
- # 1)'BayesOptDesign' 2)'BayesActDesign' 3)'VarOptDesign' 4)'alphabetic'
+ # 1)'BayesOptDesign' 2)'BayesActDesign' 3)'VarOptDesign' 4)'alphabetic'
# 5)'Space-filling'
MetaModelOpts.ExpDesign.ExploitMethod = 'VarOptDesign'
-
+
# BayesOptDesign -> when data is available
# 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
# 3)APP (A-Posterior-percision) # ['DKL', 'BME', 'infEntropy']
# MetaModelOpts.ExpDesign.UtilityFunction = 'DKL'
-
+
# VarBasedOptDesign -> when data is not available
# Only with Vornoi >>> 1)Entropy 2)EIGF, 3)LOOCV
- MetaModelOpts.ExpDesign.UtilityFunction = 'Entropy' #['EIGF', 'Entropy', 'LOOCV']
-
+ # or a combination as a list
+ MetaModelOpts.ExpDesign.UtilityFunction = 'Entropy'
+
# alphabetic
# 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
- # 3)K-Opt (K-Optimality)
- # MetaModelOpts.ExpDesign.UtilityFunction = 'D-Opt' #['D-Opt', 'A-Opt', 'K-Opt']
-
- # For colculation of validation error for SeqDesign
- MetaModelOpts.validSamples = np.load("data/Prior_"+str(ndim)+".npy")
- MetaModelOpts.validModelRuns = {'Z':np.load("data/origModelOutput_"+str(ndim)+".npy")}
- MetaModelOpts.validLikelihoods = np.load("data/validLikelihoods_"+str(ndim)+".npy")
-
+ # 3)K-Opt (K-Optimality) or a combination as a list
+ # MetaModelOpts.ExpDesign.UtilityFunction = 'D-Opt'
+
+ # For colculation of validation error for SeqDesign
+ prior = np.load(f"data/Prior_{ndim}.npy")
+ prior_outputs = np.load(f"data/origModelOutput_{ndim}.npy")
+ likelihood = np.load(f"data/validLikelihoods_{ndim}.npy")
+ MetaModelOpts.validSamples = prior
+ MetaModelOpts.validModelRuns = {'Z': prior_outputs}
+ MetaModelOpts.validLikelihoods = likelihood
+
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
- # Adaptive sparse arbitrary polynomial chaos expansion
+ # Train the meta model
PCEModel = MetaModelOpts.create_metamodel(Model)
-
- # # Save PCE models
- # with open('PCEModel_'+Model.Name+'.pkl', 'wb') as output:
- # joblib.dump(PCEModel, output, 2)
-
+
+ # Save PCE models
+ with open(f'PCEModel_{Model.name}.pkl', 'wb') as output:
+ joblib.dump(PCEModel, output, 2)
+
# Load PCEModel
- # with open('PCEModel_'+Model.Name+'.pkl', 'rb') as input:
+ # with open(f'PCEModel_{Model.name}.pkl', 'rb') as input:
# PCEModel = joblib.load(input)
- #=====================================================
- #========= POST PROCESSING OF METAMODELS ===========
- #=====================================================
+ # =====================================================
+ # ========= POST PROCESSING OF METAMODELS ===========
+ # =====================================================
PostPCE = PostProcessing(PCEModel)
-
+
# Plot to check validation visually.
PostPCE.validMetamodel(nValidSamples=3)
# Compute the moments and compare with the Monte-Carlo reference
if MetaModelOpts.metaModel != 'GPE':
PostPCE.plotMoments()
-
+
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential':
refBME_KLD = np.load("data/refBME_KLD_"+str(ndim)+".npy")
PostPCE.seqDesignDiagnosticPlots(refBME_KLD)
-
+
# Plot the sobol indices
if MetaModelOpts.metaModel != 'GPE':
sobol_cell, total_sobol = PostPCE.sobolIndicesPCE()
-
+
# Compute and print RMSE error
PostPCE.accuracyCheckMetaModel(nSamples=200)
-
- #=====================================================
- #======== Bayesian inference with Emulator ==========
- #=====================================================
+ sys.exit()
+ # =====================================================
+ # ======== Bayesian inference with Emulator ==========
+ # =====================================================
BayesOpts = BayesInference(PCEModel)
BayesOpts.emulator = True
BayesOpts.PlotPostPred = True
-
+
# BayesOpts.selectedIndices = [0, 3, 5, 7, 9]
# BME Bootstrap
# BayesOpts.Bootstrap = True
# BayesOpts.BootstrapItrNr = 500
# BayesOpts.BootstrapNoise = 100
-
+
# Bayesian cross validation
# BayesOpts.BayesLOOCV = True
-
+
# Select the inference method
BayesOpts.SamplingMethod = 'MCMC'
BayesOpts.MCMCnwalkers = 30
- # BayesOpts.MCMCinitSamples = np.zeros((1,10)) + 1e-3 * np.random.randn(200, 10)
+ # BayesOpts.MCMCinitSamples = np.zeros((1,10))+1e-3*np.random.randn(200, 10)
# BayesOpts.MCMCnSteps = 1000
import emcee
BayesOpts.MCMCmoves = emcee.moves.KDEMove()
# BayesOpts.MultiProcessMCMC = True
-
+
# ----- Define the discrepancy model -------
- # obsData = pd.DataFrame(Model.Observations, columns=Model.Output.Names)
BayesOpts.MeasurementError = obsData
-
+
# # -- (Option B) --
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
DiscrepancyOpts.Parameters = obsData**2
# DiscrepancyOpts.Parameters = {'Z':np.ones(10)}
BayesOpts.Discrepancy = DiscrepancyOpts
-
+
# -- (Option C) --
# DiscOutputOpts = Input()
# # # OutputName = 'Z'
@@ -284,7 +287,7 @@ if __name__ == "__main__":
# DiscOutputOpts.Marginals[0].Parameters = [0, 10]
# BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
# 'infer': Discrepancy(DiscOutputOpts)}
-
+
# BayesOpts.BiasInputs = {'Z':np.arange(0, 10, 1.).reshape(-1,1) / 9}
# DiscOutputOpts = Input()
# # OutputName = 'lambda'
@@ -292,78 +295,59 @@ if __name__ == "__main__":
# DiscOutputOpts.Marginals[0].Name = '$\lambda$'
# DiscOutputOpts.Marginals[0].DistType = 'unif'
# DiscOutputOpts.Marginals[0].Parameters = [0, 1]
-
+
# # OutputName = 'sigma_f'
# DiscOutputOpts.addMarginals()
# DiscOutputOpts.Marginals[1].Name = '$\sigma_f$'
# DiscOutputOpts.Marginals[1].DistType = 'unif'
# DiscOutputOpts.Marginals[1].Parameters = [0, 1e-4]
# BayesOpts.Discrepancy = Discrepancy(DiscOutputOpts)
- # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
+ # BayesOpts.Discrepancy = {'known': DiscrepancyOpts,
# 'infer': Discrepancy(DiscOutputOpts)}
-
-
+ # Start the calibration/inference
Bayes_PCE = BayesOpts.create_Inference()
-
+
# Save class objects
- with open('Bayes_'+Model.Name+'.pkl', 'wb') as output:
+ with open(f'Bayes_{Model.name}.pkl', 'wb') as output:
joblib.dump(Bayes_PCE, output, 2)
sys.exit('STOP!!')
- #=====================================================
- #====== Bayesian inference with Forward Model ========
- #=====================================================
+ # =====================================================
+ # ====== Bayesian inference with Forward Model ========
+ # =====================================================
if ndim < 7:
PCEModel.ModelObj.Name = 'OrigAnalyticFunc'
-
+
BayesOptsModel = BayesInference(PCEModel)
BayesOptsModel.emulator = False
- BayesOptsModel.NrofSamples = 10000#00
-
- # Evaluation of forward model predictions for the samples used for PCEModel
+ BayesOptsModel.NrofSamples = 10000
+
+ # Evaluation of forward model predictions for the surrogate posterior
# BayesOptsModel.Samples = Bayes_PCE.Samples
# BayesOptsModel.SamplingMethod = 'MCMC'
-
+
# Bootstrap for BME calulations
BayesOptsModel.Bootstrap = True
BayesOptsModel.BootstrapItrNr = 1
BayesOptsModel.BootstrapNoise = 0.05
-
+
BayesOptsModel.PlotPostDist = True
BayesOptsModel.PlotPostPred = False
-
-
+
# ----- Define the discrepancy model -------
BayesOptsModel.MeasurementError = obsData
# (Option B)
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = obsData **2
-
- BayesOptsModel.Discrepancy = DiscrepancyOpts
-
+ DiscrepancyOpts.Parameters = obsData ** 2
+ BayesOptsModel.Discrepancy = DiscrepancyOpts
+
+ # Start the infernce
Bayes_Model = BayesOptsModel.create_Inference()
-
- #=====================================================
- #============== Save class objects =================
- #=====================================================
- with open('AnalyticFunc_Results.pkl', 'wb') as output:
- pickle.dump(PCEModel, output, pickle.HIGHEST_PROTOCOL)
-
- pickle.dump(PostPCE, output, pickle.HIGHEST_PROTOCOL)
-
- pickle.dump(Bayes_PCE, output, pickle.HIGHEST_PROTOCOL)
-
- if ndim < 7:
- pickle.dump(Bayes_Model, output, pickle.HIGHEST_PROTOCOL)
-
-# del PCEModel
-# del PostPCE
-# del Bayes
-
- # Load the objects
-# with open('AnalyticFunc_Results.pkl', 'rb') as input:
-# PCEModel = pickle.load(input)
-# PostPCE = pickle.load(input)
-# Bayes_PCE = pickle.load(input)
-# Bayes_Model = pickle.load(input)
\ No newline at end of file
+
+ # =====================================================
+ # ============== Save class objects =================
+ # =====================================================
+ if ndim < 7:
+ with open('AnalyticFunc_Results.pkl', 'wb') as output:
+ joblib.dump(Bayes_Model, output, 2)