diff --git a/.gitignore b/.gitignore
index 367eb9c192fc14bde9fd6d749255f0b91301d2ec..194332a8e0168889d72626e3dd037f25d9c2624a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
.Rhistory
.spyproject/
+*.hdf5
BayesValidRox/.spyderproject/
# BayesInference
diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 4897cb139aa9e5712c6762a2d296a5c8135d3cdb..7b2592354c1e1bae4d6b0b2183e2c358b84a8c9a 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -679,7 +679,7 @@ class BayesInference:
else:
plotname = '/Poisterior_Dist_'+Model.Name
- figPosterior.savefig('./'+OutputDir+ plotname + '.svg',
+ figPosterior.savefig('./'+OutputDir+ plotname + '.pdf',
bbox_inches='tight')
# -------- Plot logBME dist --------
@@ -706,7 +706,7 @@ class BayesInference:
else:
plotname = '/BME_hist_'+Model.Name
- plt.savefig('./'+OutputDir+plotname+'.svg', bbox_inches='tight')
+ plt.savefig('./'+OutputDir+plotname+'.pdf', bbox_inches='tight')
plt.show()
plt.close()
@@ -846,12 +846,12 @@ class BayesInference:
- # Save figure in svg format
+ # Save figure in pdf format
if self.emulator:
plotname = '/Post_Prior_Perd_'+Model.Name+'_emulator'
else:
plotname = '/Post_Prior_Perd_'+Model.Name
- fig.savefig('./'+OutputDir+plotname+ '_' +OutputName +'.svg', bbox_inches='tight')
+ fig.savefig('./'+OutputDir+plotname+ '_' +OutputName +'.pdf', bbox_inches='tight')
return self
diff --git a/BayesValidRox/BayesInference/MCMC.py b/BayesValidRox/BayesInference/MCMC.py
index a3bf078420cff0cc060e163b82440122b4d5792a..fd8c1ec361190ba5e821c56d92eadd3ad5f8c894 100755
--- a/BayesValidRox/BayesInference/MCMC.py
+++ b/BayesValidRox/BayesInference/MCMC.py
@@ -11,6 +11,7 @@ import numpy as np
import emcee
import pandas as pd
import matplotlib.pyplot as plt
+from matplotlib.backends.backend_pdf import PdfPages
from multiprocessing import Pool
import scipy.stats as st
from scipy.linalg import cholesky as chol
@@ -107,7 +108,10 @@ class MCMC():
params_range = PCEModel.ExpDesign.BoundTuples
OutputDir = (r'Outputs_Bayes_' + Model.Name + '_' + self.BayesOpts.Name)
if not os.path.exists(OutputDir): os.makedirs(OutputDir)
-
+
+ # create a PdfPages object
+ pdf = PdfPages(OutputDir+'/traceplots.pdf')
+ fig = plt.figure()
for parIdx in range(ndim):
# Set up the axes with gridspec
fig = plt.figure()
@@ -128,7 +132,13 @@ class MCMC():
main_ax.set_title('traceplot for ' + parNames[parIdx])
main_ax.set_xlabel('step number')
- fig.savefig(OutputDir+'/traceplot_par_'+str(parIdx+1)+'.svg', bbox_inches='tight')
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
+
+ pdf.close()
# logml_dict = self.Marginal_llk_emcee(sampler, self.nburn, logp=None, maxiter=5000)
# print('\nThe Bridge Sampling Estimation is %.5f.'%(logml_dict['logml']))
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 21e2bb6e4719a83ac85037f470e02e6ed8516aca..e7356e3163249c8c040d5a3d6280574925100661 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -10,6 +10,7 @@ import numpy as np
import math
import os
from itertools import combinations, cycle
+from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.rcParams.update({'font.size': 24})
@@ -40,25 +41,24 @@ class PostProcessing:
self.ModelOutputs = []
self.PCEOutputs = []
self.PCEOutputs_std = []
- self.plotFlag = False
self.sobol_cell = {}
self.total_sobol = {}
self.Likelihoods = []
self.Name = 'calib'
#--------------------------------------------------------------------------------------------------------
- def PCEMoments(self, aPCE):
+ def PCEMoments(self):
PCEModel = self.PCEModel
Model = PCEModel.ModelObj
- for Outkey, ValuesDict in aPCE.CoeffsDict.items():
+ for Outkey, ValuesDict in PCEModel.CoeffsDict.items():
PCEMean = np.zeros((len(ValuesDict)))
PCEVar = np.zeros((len(ValuesDict)))
for Inkey, InIdxValues in ValuesDict.items():
idx = int(Inkey.split('_')[1]) - 1
- coeffs = aPCE.CoeffsDict[Outkey][Inkey]
+ coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
# Mean = c_0
PCEMean[idx] = coeffs[0]
@@ -76,51 +76,69 @@ class PostProcessing:
self.Reference = Model.read_MCReference()
#--------------------------------------------------------------------------------------------------------
- def MomentLineplot(self, PCE_Means, PCE_Stds , x_label, SaveFig = False):
+ def plotMoments(self, x_label='Time [s]', SaveFig=True):
+
+ # Compute the moments with the PCEModel object
+ self.PCEMoments()
+
+ # Open a pdf for the plots
+ if SaveFig:
+ newpath = (r'Outputs_PostProcessing')
+ if not os.path.exists(newpath): os.makedirs(newpath)
+
+ # create a PdfPages object
+ pdf = PdfPages('./'+newpath+'/Mean_Std_PCE.pdf')
# Create the plot object
- NrofOutputs = len(PCE_Means)
- Keys = list(PCE_Means.keys())
+ NrofOutputs = len(self.PCEMeans)
+ Keys = list(self.PCEMeans.keys())
fig, ax = plt.subplots(nrows=NrofOutputs, ncols=2)
if NrofOutputs == 1:
ax = ax[:,None].T
- x_data = self.Reference[x_label]
+ # Set the x values
+ x_values = self.PCEModel.ExpDesign.Y['x_values']
+ metaModel = self.PCEModel.metaModel
# Plot the best fit line, set the linewidth (lw), color and
# transparency (alpha) of the line
for idx in range(NrofOutputs):
- mean_data = PCE_Means[Keys[idx]]
- std_data = PCE_Stds[Keys[idx]]
+ mean_data = self.PCEMeans[Keys[idx]]
+ std_data = self.PCEStd[Keys[idx]]
- ax[idx, 0].plot(x_data, mean_data, lw = 2, color = '#539caf', alpha = 1, label = 'BaPCE')
- ax[idx, 0].plot(x_data, self.Reference['mean'], lw = 2, color = 'r', alpha = 1, label = 'Ref.')
- ax[idx, 1].plot(x_data, std_data, lw = 2, color = '#539caf', alpha = 1, label = 'BaPCE')
- ax[idx, 1].plot(x_data, self.Reference['std'], lw = 2, color = 'r', alpha = 1, label = 'Ref.')
+ ax[idx,0].plot(x_values, mean_data, lw=3, color='#539caf', marker='x', label=metaModel)
+ ax[idx,1].plot(x_values, std_data, lw=3, color='#539caf', marker='x', label=metaModel)
+
+ try:
+ ax[idx,0].plot(x_values, self.Reference['mean'], lw=3, marker='x', color='r', label='Ref.')
+ ax[idx,1].plot(x_values, self.Reference['std'], lw=3, marker='x', color='r', label='Ref.')
+ except:
+ pass
# Label the axes and provide a title
- ax[idx, 0].set_xlabel(x_label)
- ax[idx, 1].set_xlabel(x_label)
- ax[idx, 0].set_ylabel('Mean of '+ Keys[idx])
- ax[idx, 1].set_ylabel('Std. of '+ Keys[idx])
+ ax[idx,0].set_xlabel(x_label)
+ ax[idx,1].set_xlabel(x_label)
+ ax[idx,0].set_ylabel(Keys[idx])
+ ax[idx,1].set_ylabel(Keys[idx])
+
+ # Provide a title
+ ax[0,0].set_title('Mean of '+ Keys[idx])
+ ax[0,1].set_title('Std of '+ Keys[idx])
+
+ ax[0,0].legend(loc='best')
+ ax[0,1].legend(loc='best')
+
+ plt.tight_layout()
+
+ if SaveFig:
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
+
+ pdf.close()
- # Provide a title
- ax[0, 0].set_title("Mean")
- ax[0, 1].set_title("Std")
-
- ax[0, 0].legend(loc='best')
- ax[0, 1].legend(loc='best')
-
- plt.tight_layout()
-
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = (r'Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- plt.savefig('./'+newpath+'/Mean_Std_PCE.svg') # save the figure to file
- plt.close(fig)
-
-
#--------------------------------------------------------------------------------------------------------
def get_Sample(self):
@@ -135,6 +153,8 @@ class PostProcessing:
#--------------------------------------------------------------------------------------------------------
def eval_PCEmodel_3D(self, SaveFig=True):
+ self.NrofSamples = 1000
+
PCEModel = self.PCEModel
Model = self.PCEModel.ModelObj
NrofSamples = self.NrofSamples
@@ -180,7 +200,7 @@ class PostProcessing:
PCEOutputs_std[:, idx] = y_std
# Model evaluation
- ModelOutputsDict, _ = Model.Run_Model_Parallel(SampleMesh)
+ ModelOutputsDict, _ = Model.Run_Model_Parallel(SampleMesh, keyString='Valid3D')
ModelOutputs[:,idx]= ModelOutputsDict[Outkey].T
@@ -206,7 +226,7 @@ class PostProcessing:
if not os.path.exists(newpath): os.makedirs(newpath)
# 3D-plot PCEModel
- fig_PCE.savefig('./'+newpath+'/3DPlot_PCEModel.svg', format="svg",
+ fig_PCE.savefig('./'+newpath+'/3DPlot_PCEModel.pdf', format="pdf",
bbox_inches='tight') # save the figure to file
plt.close(fig_PCE)
@@ -227,7 +247,7 @@ class PostProcessing:
if SaveFig is True:
# ---------------- Saving the figure and text files -----------------------
# 3D-plot Model
- fig_Model.savefig('./'+newpath+'/3DPlot_Model.svg', format="svg",
+ fig_Model.savefig('./'+newpath+'/3DPlot_Model.pdf', format="pdf",
bbox_inches='tight') # save the figure to file
plt.close(fig_Model)
@@ -236,7 +256,7 @@ class PostProcessing:
#--------------------------------------------------------------------------------------------------------
- def eval_Model(self, Samples=None):
+ def eval_Model(self, Samples=None, keyString='Valid'):
"""
Evaluate Forward Model
@@ -251,7 +271,7 @@ class PostProcessing:
Samples = Samples
self.NrofSamples = len(Samples)
- ModelOutputs, CollocationPoints = Model.Run_Model_Parallel(Samples, keyString='Valid')
+ ModelOutputs, CollocationPoints = Model.Run_Model_Parallel(Samples, keyString)
if index is not None:
self.ModelOutputs = {key: ModelOutputs[key][:index] for key in [list(ModelOutputs.keys())[0]]}
@@ -264,32 +284,33 @@ class PostProcessing:
#--------------------------------------------------------------------------------------------------------
- def ValidMetamodel(self, x_values=[], x_axis="x [m]"):
+ def validMetamodel(self, nValidSamples=1, x_axis="Time [s]"):
"""
Evaluate the meta model and the PCEModel
"""
metaModel = self.PCEModel
+ self.NrofSamples = nValidSamples
Samples = self.get_Sample()
+ x_values = self.PCEModel.ExpDesign.Y['x_values']
- self.eval_Model(Samples)
+ self.eval_Model(Samples, keyString='valid')
self.PCEOutputs, self.PCEOutputs_std = metaModel.eval_metamodel(samples=Samples,name=self.Name)
- if self.plotFlag == True:
- try:
- key = list(self.ModelOutputs.keys())[1]
- except:
- key = list(self.ModelOutputs.keys())[0]
-
- NrofOutputs = self.ModelOutputs[key].shape[1]
+ try:
+ key = list(self.ModelOutputs.keys())[1]
+ except:
+ key = list(self.ModelOutputs.keys())[0]
- if NrofOutputs == 1:
- self.plotValidation()
- else:
- self.plotValidationMulti(x_values=x_values, x_axis=x_axis)
+ NrofOutputs = self.ModelOutputs[key].shape[1]
+
+ if NrofOutputs == 1:
+ self.plotValidation()
+ else:
+ self.plotValidationMulti(x_values=x_values, x_axis=x_axis)
#--------------------------------------------------------------------------------------------------------
- def relErrorPCEModel(self, nSamples=None, Samples=None, validOutputsDict=None):
+ def accuracyCheckMetaModel(self, nSamples=None, Samples=None, validOutputsDict=None):
"""
Evaluate the relative error of the PCEModel
"""
@@ -302,7 +323,7 @@ class PostProcessing:
Samples = self.get_Sample() if Samples is None else Samples
# Run the original model with the generated samples
- ModelOutputs = self.eval_Model(Samples) if validOutputsDict is None else validOutputsDict
+ ModelOutputs = self.eval_Model(Samples, keyString='validSet') if validOutputsDict is None else validOutputsDict
# Run the PCE model with the generated samples
PCEOutputs, PCEOutputs_std = metaModel.eval_metamodel(samples=Samples)
@@ -323,8 +344,7 @@ class PostProcessing:
1) https://towardsdatascience.com/how-do-you-check-the-quality-of-your-regression-model-in-python-fa61759ff685
"""
- aPCE = self.PCEModel
-
+ PCEModel = self.PCEModel
# get the samples
@@ -334,9 +354,21 @@ class PostProcessing:
Y_PC_Val = self.PCEOutputs
Y_Val = self.ModelOutputs
+ # Open a pdf for the plots
+ if SaveFig:
+ newpath = (r'Outputs_PostProcessing')
+ if not os.path.exists(newpath): os.makedirs(newpath)
+
+ # create a PdfPages object
+ pdf1 = PdfPages('./'+newpath+'/Model_vs_PCEModel.pdf')
+ pdf2 = PdfPages('./'+newpath+'/Residuals_vs_PredVariables.pdf')
+ pdf3 = PdfPages('./'+newpath+'/Fitted_vs_Residuals.pdf')
+ pdf4 = PdfPages('./'+newpath+'/Hist_NormResiduals.pdf')
+ pdf5 = PdfPages('./'+newpath+'/QQPlot_NormResiduals.pdf')
+
+ fig = plt.figure()
# Fit the data(train the model)
for key in Y_PC_Val.keys():
- fig = plt.figure()
Y_PC_Val_ = Y_PC_Val[key]
Y_Val_ = Y_Val[key]
@@ -353,7 +385,7 @@ class PostProcessing:
# Calculate the adjusted R_squared and RMSE
# the total number of explanatory variables in the model (not including the constant term)
- length_list = [len(value) for Key, value in aPCE.BasisDict[key].items()]
+ length_list = [len(value) for Key, value in PCEModel.BasisDict[key].items()]
NofPredictors = min(length_list)
TotalSampleSize = X_Val.shape[0] #sample size
@@ -369,18 +401,16 @@ class PostProcessing:
plt.grid()
plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = ('Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig.savefig('./'+newpath+'/Model_vs_PCEModel_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
- plt.close(fig)
+ if SaveFig:
+ # save the current figure
+ pdf1.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
# ------ Residuals vs. predicting variables ------
# Check the assumptions of linearity and independence
- fig1 = plt.figure()
- plt.title("Residuals vs. predicting variables" ,fontsize=16)
+ plt.title(key+": Residuals vs. predicting variables" ,fontsize=16)
residuals = Y_Val_ - Y_PC_Val_
plt.scatter(x=Y_Val_,y=residuals,color='blue',edgecolor='k')
plt.grid(True)
@@ -390,19 +420,18 @@ class PostProcessing:
plt.xlabel(key,fontsize=14)
plt.ylabel('Residuals',fontsize=14)
plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = ('Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig1.savefig('./'+newpath+'/Residuals_vs_PredVariables_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
- plt.close(fig1)
+ if SaveFig:
+ # save the current figure
+ pdf2.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
# ------ Fitted vs. residuals ------
# Check the assumptions of linearity and independence
- fig2 = plt.figure()
- plt.title("Residuals vs. predicting variables" ,fontsize=16)
+ # fig2 = plt.figure()
+ plt.title(key+": Residuals vs. predicting variables" ,fontsize=16)
residuals = Y_Val_ - Y_PC_Val_
plt.scatter(x=Y_PC_Val_,y=residuals,color='blue',edgecolor='k')
plt.grid(True)
@@ -412,21 +441,21 @@ class PostProcessing:
plt.xlabel(key,fontsize=14)
plt.ylabel('Residuals',fontsize=14)
plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = ('Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig2.savefig('./'+newpath+'/Fitted_vs_Residuals_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
- plt.close(fig2)
-
+
+ if SaveFig:
+ # save the current figure
+ pdf3.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
+
# ------ Histogram of normalized residuals ------
- fig3 = plt.figure()
+ # fig3 = plt.figure()
resid_pearson = residuals / (max(residuals)-min(residuals))
plt.hist(resid_pearson,bins=20,edgecolor='k')
plt.ylabel('Count')
plt.xlabel('Normalized residuals')
- plt.title("Histogram of normalized residuals",fontsize=18)
+ plt.title(key+": Histogram of normalized residuals",fontsize=18)
# Normality (Shapiro-Wilk) test of the residuals
ax = plt.gca()
@@ -445,13 +474,13 @@ class PostProcessing:
ax.add_artist(at)
plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = ('Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig3.savefig('./'+newpath+'/Hist_NormResiduals_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
- plt.close(fig3)
+
+ if SaveFig:
+ # save the current figure
+ pdf4.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
# ------ Q-Q plot of the normalized residuals ------
from statsmodels.graphics.gofplots import qqplot
@@ -461,25 +490,36 @@ class PostProcessing:
plt.yticks()
plt.xlabel("Theoretical quantiles")
plt.ylabel("Sample quantiles")
- plt.title("Q-Q plot of normalized residuals",fontsize=18)
+ plt.title(key+": Q-Q plot of normalized residuals",fontsize=18)
plt.grid(True)
plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = ('Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig4.savefig('./'+newpath+'/QQPlot_NormResiduals_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
- plt.close(fig4)
-
+ if SaveFig:
+ # save the current figure
+ pdf5.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
+
+ # Close the pdfs
+ pdf1.close()
+ pdf2.close()
+ pdf3.close()
+ pdf4.close()
+ pdf5.close()
#--------------------------------------------------------------------------------------------------------
def plotValidationMulti(self, x_values=[], x_axis="x [m]", SaveFig=True):
Model = self.PCEModel.ModelObj
-
+ if SaveFig:
+ newpath = (r'Outputs_PostProcessing')
+ if not os.path.exists(newpath): os.makedirs(newpath)
+
+ # create a PdfPages object
+ pdf = PdfPages('./'+newpath+'/Model_vs_PCEModel.pdf')
+
# List of markers and colors
color = cycle((['b', 'g', 'r', 'y', 'k']))
marker = cycle(('x', 'd', '+', 'o', '*'))
@@ -496,10 +536,11 @@ class PostProcessing:
x_values = Y_Val[OutNames[0]]
except:
x_values = x_values
-
+
+ fig = plt.figure(figsize=(24, 16))
+
# Plot the model vs PCE model
for keyIdx, key in enumerate(OutNames[1:]):
- fig = plt.figure(figsize=(24, 16))
Y_PC_Val_ = Y_PC_Val[key]
Y_PC_Val_std_ = Y_PC_Val_std[key]
@@ -512,10 +553,17 @@ class PostProcessing:
plt.plot(x_values, Y_Val_[idx,:], color=Color, marker=Marker, label='$Y_{%s}^{M}$'%(idx+1))
plt.plot(x_values, Y_PC_Val_[idx,:], color=Color, marker=Marker, linestyle='--', label='$Y_{%s}^{PCE}$'%(idx+1))
- #plt.errorbar(x_values, Y_PC_Val_[idx,:], 2*Y_PC_Val_std_[idx,:], color=Color, marker=Marker, linestyle='--')
plt.fill_between(x_values, Y_PC_Val_[idx,:]-1.96*Y_PC_Val_std_[idx,:],
Y_PC_Val_[idx,:]+1.96*Y_PC_Val_std_[idx,:], color=Color,alpha=0.15)
-
+
+ # Calculate the RMSE
+ RMSE = mean_squared_error(Y_PC_Val_, Y_Val_, squared=False)
+ R2 = r2_score(Y_Val_[idx,:].reshape(-1,1), Y_Val_[idx,:].reshape(-1,1))
+
+ plt.annotate('RMSE = '+ str(round(RMSE, 3)) + '\n' + r'$R^2$ = '+ str(round(R2, 3)),
+ xy=(0.2, 0.75), xycoords='axes fraction')
+
+
plt.ylabel(key)
plt.xlabel(x_axis)
plt.legend(loc='best')
@@ -523,14 +571,14 @@ class PostProcessing:
#plt.show()
- if SaveFig is True:
- # ---------------- Saving the figure and text files -----------------------
- newpath = (r'Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- plt.savefig('./'+newpath+'/Model_vs_PCEModel_'+key+'.svg',
- bbox_inches='tight') # save the figure to file
-
- plt.close(fig)
+ if SaveFig:
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
+
+ pdf.close()
# Cleanup
#Zip the subdirectories
@@ -543,7 +591,7 @@ class PostProcessing:
#--------------------------------------------------------------------------------------------------------
- def seqDesignDiagnosticPlots(self, NrofInitSamples, refBME_KLD=None, SaveFig=False):
+ def seqDesignDiagnosticPlots(self, refBME_KLD=None, SaveFig=True):
"""
Plot the Kullback-Leibler divergence in the sequential design.
@@ -552,6 +600,13 @@ class PostProcessing:
NrofInitSamples = PCEModel.ExpDesign.initNrSamples
totalNSamples = PCEModel.ExpDesign.X.shape[0]
+ if SaveFig:
+ newpath = (r'Outputs_PostProcessing')
+ if not os.path.exists(newpath): os.makedirs(newpath)
+
+ # create a PdfPages object
+ pdf = PdfPages('./'+newpath+'/seqPCEModelDiagnostics.pdf')
+
plotList = ['Modified LOO error', 'Validation error', 'KLD', 'BME',
'RMSEMean', 'RMSEStd', 'Hellinger distance']
seqList = [PCEModel.SeqModifiedLOO,PCEModel.seqValidError, PCEModel.SeqKLD,
@@ -561,6 +616,7 @@ class PostProcessing:
markers = ('x', 'o', 'd', '*', '+')
colors = ('k', 'darkgreen', 'b', 'navy', 'darkred')
+ fig, ax = plt.subplots()
# Plot the evolution of the diagnostic criteria of the Sequential Experimental Design.
for plotidx, plot in enumerate(plotList):
@@ -614,7 +670,7 @@ class PostProcessing:
step2 = 1 if PCEModel.ExpDesign.NrofNewSample!=1 else 5
edge_color = ['red', 'blue', 'green']
fill_color = ['tan', 'cyan', 'lightgreen']
- fig, ax = plt.subplots()
+
plotLabel = plot
# Plot for different Utility Functions
for idx, util in enumerate(UtilFuncs):
@@ -657,24 +713,24 @@ class PostProcessing:
plt.autoscale(True)
ax.set_xlabel('Size of the experimental design')
ax.set_ylabel(plotLabel)
+ plt.title(plot)
- if SaveFig is True:
- # ---------------- Saving the figure -----------------------
- newpath = (r'Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig.savefig('./'+newpath+'/seqPCEModel_'+plot+'.svg',
- bbox_inches='tight') # save the figure to file
+ if SaveFig:
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
- # ---------------- Saving arrays into files -----------------------
+ # Destroy the current plot
+ plt.clf()
+
+ # Save arrays into files
#np.savetxt('./'+newpath+'/Seq'+plot+'.txt', sortedSeqOpt)
f = open('./'+newpath+'/Seq'+plot+'.txt',"w")
f.write( str(sortedSeqOpt) )
f.close()
- plt.close(fig)
else:
- fig,ax = plt.subplots()
+ # fig,ax = plt.subplots()
for idx, name in enumerate(nameUtil):
SeqValues = SeqDict[name]
@@ -730,26 +786,26 @@ class PostProcessing:
plt.xlabel('Number of runs')
plt.ylabel(plotLabel)
+ plt.title(plot)
plt.legend()
- if SaveFig is True:
- # ---------------- Saving the figure -----------------------
- newpath = (r'Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- plt.savefig('./'+newpath+'/seqPCEModel_'+plot+'.svg',
- bbox_inches='tight') # save the figure to file
+ if SaveFig:
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
+ # Destroy the current plot
+ plt.clf()
+
# ---------------- Saving arrays into files -----------------------
np.save('./'+newpath+'/Seq'+plot+'.npy', SeqValues)
- plt.close(fig)
-
-
+ # Close the pdf
+ pdf.close()
return
#--------------------------------------------------------------------------------------------------------
- def SobolIndicesPCE(self, x_values, x_axis, barWidth=0.1, colormap='gray', PlotType='bar', SaveFig = False):
+ def sobolIndicesPCE(self, x_axis='Time [s]', PlotType=None, SaveFig=True):
"""
Gives information about the importance of parameters
More Info: Eq. 27 @ Global sensitivity analysis: A flexible and efficient framework with an example from stochastic hydrogeology
@@ -765,6 +821,17 @@ class PostProcessing:
PCEModel = self.PCEModel
NofPa = PCEModel.NofPa
max_order = PCEModel.MaxPceDegree
+ x_values = self.PCEModel.ExpDesign.Y['x_values']
+
+ if SaveFig:
+ newpath = (r'Outputs_PostProcessing')
+ if not os.path.exists(newpath): os.makedirs(newpath)
+
+ # create a PdfPages object
+ pdf = PdfPages('./'+newpath+'/Sobol_indices.pdf')
+
+ fig = plt.figure()
+
for Output in PCEModel.ModelObj.Output.Names:
@@ -912,26 +979,29 @@ class PostProcessing:
self.total_sobol[Output] = total_sobol_array
# ---------------- Plot -----------------------
- fig = plt.figure()
parNames = [PCEModel.Inputs.Marginals[i].Name for i in range(NofPa)]
- for i, sobolIndices in enumerate(self.total_sobol[Output]):
- plt.plot(x_values, sobolIndices, label=parNames[i], lw=2.5)
-
- plt.ylabel('Total Sobol indices, $S^T$')
- plt.xlabel(x_axis)
+ if PlotType =='bar':
+ ax = fig.add_axes([0,0,1,1])
+ ax.bar(parNames,self.total_sobol[Output].reshape(1,-1)[0])
+ else:
+ for i, sobolIndices in enumerate(self.total_sobol[Output]):
+ plt.plot(x_values, sobolIndices, label=parNames[i], marker='x', lw=2.5)
+
+ plt.ylabel('Total Sobol indices, $S^T$')
+ plt.xlabel(x_axis)
+
plt.title('Sensitivity analysis of '+ Output)
plt.legend(loc='best')
- if SaveFig is True:
- # ---------------- Saving the figure -----------------------
- newpath = (r'Outputs_PostProcessing')
- if not os.path.exists(newpath): os.makedirs(newpath)
- fig.savefig('./'+newpath+'/Sobol_indices_'+Output+'.svg', bbox_inches='tight') # save the figure to file
+ if SaveFig:
+ # save the current figure
+ pdf.savefig(fig, bbox_inches='tight')
- try:
- plt.close(fig)
- except: pass
-
+ # Destroy the current plot
+ plt.clf()
+
+ pdf.close()
+
return self.sobol_cell, self.total_sobol
\ No newline at end of file
diff --git a/BayesValidRox/PostProcessing/adaptPlot.py b/BayesValidRox/PostProcessing/adaptPlot.py
index 080462bacaa47e755fd3cfa7b70c41063211bc89..5b9151b9174d583a034f58b416123aedb64ad469 100755
--- a/BayesValidRox/PostProcessing/adaptPlot.py
+++ b/BayesValidRox/PostProcessing/adaptPlot.py
@@ -81,11 +81,13 @@ def adaptPlot(PCEModel, Y_Val, Y_PC_Val, Y_PC_Val_std, x_values=[], plotED=False
RMSE = mean_squared_error(Y_PC_Val_, Y_Val_, squared=False)
R2 = r2_score(Y_PC_Val_.reshape(-1,1), Y_Val_.reshape(-1,1))
- plt.annotate('RMSE = '+ str(round(RMSE, 3)) + '\n' + r'$R^2$ = '+ str(round(R2, 3)), xy=(0.05, 0.85), xycoords='axes fraction')
+ plt.annotate('RMSE = '+ str(round(RMSE, 3)) + '\n' + r'$R^2$ = '+ str(round(R2, 3)),
+ xy=(0.2, 0.75), xycoords='axes fraction')
plt.ylabel(key)
plt.xlabel(x_axis)
plt.legend(loc='best')
+ plt.title(key)
plt.grid()
if SaveFig:
@@ -95,4 +97,4 @@ def adaptPlot(PCEModel, Y_Val, Y_PC_Val, Y_PC_Val_std, x_values=[], plotED=False
# Destroy the current plot
plt.clf()
- pdf.close()
\ No newline at end of file
+ pdf.close()
\ No newline at end of file
diff --git a/BayesValidRox/PyLink/FuncForwardModel.py b/BayesValidRox/PyLink/FuncForwardModel.py
index 6bd9be31d910a08d4dacbc89724ebd6c3dc7275b..d2db1b0687672019f5e03d240d8f8a52c6cf61dd 100644
--- a/BayesValidRox/PyLink/FuncForwardModel.py
+++ b/BayesValidRox/PyLink/FuncForwardModel.py
@@ -73,7 +73,7 @@ class FuncForwardModel:
grpY = f.create_group("EDY/"+var) if not hdf5_exist else f.get("EDY/"+var)
Outputs = np.asarray([item[varIdx+1] for item in group_results], dtype=np.float64)
if prevRun_No == 0:
- grpY.create_dataset("init_"+keyString, data=CollocationPoints)
+ grpY.create_dataset("init_"+keyString, data=Outputs)
else:
try:
oldEDY = np.array(f['EDY/'+var+'/adaptive_'+keyString])
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index f5b12d9e8bafce2b367ec1483ad736d336313b2b..92d72abecd5326635e1c578864a4c5b99769b9be 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -53,7 +53,6 @@ from .ExpDesigns import ExpDesigns
from .aPoly_Construction import aPoly_Construction
from .Exploration import Exploration
from .RegressionFastARD import RegressionFastARD
-from .VarLinearRegression import VarLinearRegression
class Metamodel:
@@ -546,9 +545,6 @@ class Metamodel:
compute_score=compute_score,
n_iter=1000, tol=1e-3)
- elif RegMethod.lower() == 'varbayes':
- clf_poly = VarLinearRegression(n_iter = 500, tol =1e-3)
-
elif RegMethod.lower() == 'lars':
clf_poly = linear_model.Lars(fit_intercept=True)
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
index 8ddc8cebeb81650651b606891a9c0350db30a54b..4498b6567343491e3e3597ff52d2bab1dbaf32af 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
@@ -201,32 +201,21 @@ if __name__ == "__main__":
PostPCE = PostProcessing(PCEModel)
PostPCE.Name = 'Calib'
- PostPCE.NrofSamples = 3
- PostPCE.plotFlag = True
- t = np.arange(0, 10, 1.) / 5
- PostPCE.ValidMetamodel(x_values=t, x_axis="Time [s]")
-
-# # Comparison with Monte-Carlo reference
-# # Compute the moments
-# PostPCE.PCEMoments(PCEModelCalib)
-#
-# # Plot Mean & Std for all Outputs
-# PostPCE.MomentLineplot(PCE_Means = PostPCE.PCEMeans,
-# PCE_Stds = PostPCE.PCEStd,
-# x_label = 'Time [s]', SaveFig = True)
+ PostPCE.validMetamodel(nValidSamples= 3)
+
+ # Compute the moments and compare with the Monte-Carlo reference
+# PostPCE.plotMoments()
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential':
- PostPCE.seqDesignDiagnosticPlots(NrofInitSamples, #refBME_KLD=(-19.28640,2.582278348728405)
- SaveFig=True)
+ PostPCE.seqDesignDiagnosticPlots()
# Plot the sobol indices
- sobol_cell, total_sobol = PostPCE.SobolIndicesPCE(x_values=t, x_axis="Time [s]",
- SaveFig=True)
+ sobol_cell, total_sobol = PostPCE.sobolIndicesPCE()
# Compute and print RMSE error
- PostPCE.relErrorPCEModel(nSamples=3000)
+ PostPCE.accuracyCheckMetaModel(nSamples=3000)
#=====================================================
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index fc1a2ce3653039b4814b5df29f7a3e727a7b4cb2..ddaa491e8bf4afef13064ca2da226c1e4c5afca9 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -39,7 +39,7 @@ from BayesInference.BayesInference import BayesInference, Discrepancy
if __name__ == "__main__":
# Number of parameters
- ndim = 10 # 2, 10
+ ndim = 2 # 2, 10
#=====================================================
#============= COMPUTATIONAL MODEL ================
@@ -58,7 +58,6 @@ if __name__ == "__main__":
Model.Observations['Z'] = np.repeat([2],10)
# For Checking with the MonteCarlo refrence
- #Model.MCReferenceFile = 'MCrefs_MeanStd.csv'
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")
@@ -111,7 +110,7 @@ if __name__ == "__main__":
# error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
MetaModelOpts.MinPceDegree = 1
- MetaModelOpts.MaxPceDegree = 10
+ MetaModelOpts.MaxPceDegree = 6
# q-quasi-norm 0<q<1 (default=1)
MetaModelOpts.q = 1.0 if ndim<5 else 0.75
@@ -126,7 +125,7 @@ if __name__ == "__main__":
MetaModelOpts.addExpDesign()
# One-shot (normal) or Sequential Adaptive (sequential) Design
- MetaModelOpts.ExpDesign.Method = 'sequential'
+ MetaModelOpts.ExpDesign.Method = 'normal'
NrofInitSamples = 75 #75
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
@@ -141,7 +140,7 @@ if __name__ == "__main__":
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
- MetaModelOpts.ExpDesign.MaxNSamples = 100 #150
+ MetaModelOpts.ExpDesign.MaxNSamples = 85 #150
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-16
# Defining the measurement error, if it's known a priori
@@ -208,35 +207,24 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 1
- PostPCE.plotFlag = True
- t = np.arange(0, 10, 1.) / 9
- PostPCE.ValidMetamodel(x_values=t, x_axis="Time [s]")
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=1)
- # Comparison with Monte-Carlo reference
- # Compute the moments
- PostPCE.PCEMoments(PCEModel)
-
- # Plot Mean & Std for all Outputs
- PostPCE.MomentLineplot(PCE_Means = PostPCE.PCEMeans,
- PCE_Stds = PostPCE.PCEStd,
- x_label = 'Time [s]', SaveFig = True)
-
+ # Compute the moments and compare with the Monte-Carlo reference
+ 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(NrofInitSamples, refBME_KLD
- , SaveFig=True)
+ PostPCE.seqDesignDiagnosticPlots(refBME_KLD)
# Plot the sobol indices
- sobol_cell, total_sobol = PostPCE.SobolIndicesPCE(x_values=t, x_axis="Time [s]",
- SaveFig=True)
+ sobol_cell, total_sobol = PostPCE.sobolIndicesPCE()
# Compute and print RMSE error
- PostPCE.relErrorPCEModel(nSamples=3000)
+ PostPCE.accuracyCheckMetaModel(nSamples=200)
- sys.exit('STOP!!')
+ # sys.exit('STOP!!')
#=====================================================
#======== Bayesian inference with Emulator ==========
#=====================================================
diff --git a/BayesValidRox/tests/BeamTest/Test_Beam.py b/BayesValidRox/tests/BeamTest/Test_Beam.py
index 2477b327198ab01ed5700f54f91d2c13dc6608fa..f09021897184a4c2f2356942fb304786ed865369 100644
--- a/BayesValidRox/tests/BeamTest/Test_Beam.py
+++ b/BayesValidRox/tests/BeamTest/Test_Beam.py
@@ -191,20 +191,12 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- # Compute the moments
- PostPCE.PCEMoments(PCEModel)
+ # Compute the moments and compare with the Monte-Carlo reference
+ PostPCE.plotMoments()
- # Comparison with Monte-Carlo reference
- # Plot Mean & Std for all Outputs
- PostPCE.MomentLineplot(PCE_Means = PostPCE.PCEMeans,
- PCE_Stds = PostPCE.PCEStd,
- x_label = 'x [m]', SaveFig = True)
-
- PostPCE.NrofSamples = 3
- PostPCE.plotFlag = True
-
- PostPCE.ValidMetamodel()
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=3)
#=====================================================
#============ Bayesian calibration =================
diff --git a/BayesValidRox/tests/BraninTest/Branin.py b/BayesValidRox/tests/BraninTest/Branin.py
index d831f6802f8e6d9fb1f0345de407927a53244c8c..21070e320566c50965a5ddb2fc9de1ca1a1620cf 100644
--- a/BayesValidRox/tests/BraninTest/Branin.py
+++ b/BayesValidRox/tests/BraninTest/Branin.py
@@ -27,5 +27,5 @@ def Branin(xx, *args):
Y = a*np.square(xx[:,1] - b*np.square(xx[:,0]) + c*xx[:,0] - r) + s*(1-t)*np.cos(xx[:,0]) + s
-
- return np.array(Y)[:,None]
\ No newline at end of file
+ return np.vstack((np.zeros(1),Y))
+ # return np.array(Y)[:,None]
\ No newline at end of file
diff --git a/BayesValidRox/tests/BraninTest/Branin_Test.py b/BayesValidRox/tests/BraninTest/Branin_Test.py
index 4426e5c815526ead00cfb8ad272a7115da37d874..f20d043f9323ad0b76d2c1dcc9a9521ceefa36c5 100644
--- a/BayesValidRox/tests/BraninTest/Branin_Test.py
+++ b/BayesValidRox/tests/BraninTest/Branin_Test.py
@@ -19,6 +19,7 @@ plt.rc('axes', labelsize=16)
import seaborn as sns
import os
+import pandas as pd
try:
import cPickle as pickle
@@ -27,7 +28,7 @@ except ModuleNotFoundError:
from PyLink.FuncForwardModel import FuncForwardModel
from surrogate_models.Input import Input
-from surrogate_models.surrogate_models import aPCE
+from surrogate_models.surrogate_models import Metamodel
from PostProcessing.PostProcessing import PostProcessing
from BayesInference.BayesInference import BayesInference, Discrepancy
@@ -36,7 +37,7 @@ if __name__ == "__main__":
#=====================================================
#============= COMPUTATIONAL MODEL ================
#=====================================================
- Model = FuncForwardModel('Branin')
+ Model = FuncForwardModel()
Model.Type = 'Function'
Model.pyFile = 'Branin'
@@ -44,8 +45,8 @@ if __name__ == "__main__":
Model.Output.Names = ['Z']
- Model.Observations = np.array([[45.62351]]) #[np.pi,9]
-
+ Model.Observations['Time [s]'] = np.zeros((1))
+ Model.Observations['Z'] = np.array([45.62351]) #[np.pi,9]
#=====================================================
#========= PROBABILISTIC INPUT MODEL ==============
@@ -68,7 +69,7 @@ if __name__ == "__main__":
#=====================================================
#====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
#=====================================================
- MetaModelOpts = aPCE(Inputs)
+ MetaModelOpts = Metamodel(Inputs)
# Specify the max degree to be compared by the adaptive algorithm:
# The degree with the lowest Leave-One-Out cross-validation (LOO)
@@ -78,10 +79,10 @@ if __name__ == "__main__":
# Select the sparse least-square minimization method for
# the PCE coefficients calculation:
- # 1)AaPCE: Adaptive aPCE 2)BRR: Bayesian Ridge Regression
- # 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression (Not Working)
-
- MetaModelOpts.RegMethod = 'BRR'
+ # 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
+ MetaModelOpts.RegMethod = 'FastARD'
# Print summary of the regression results
#MetaModelOpts.DisplayFlag = True
@@ -91,10 +92,15 @@ if __name__ == "__main__":
# hypercube sampling of the input model or user-defined values of X and/or Y:
MetaModelOpts.addExpDesign()
- NrofInitSamples = 3
+ # One-shot (normal) or Sequential Adaptive (sequential) Design
+ MetaModelOpts.ExpDesign.Method = 'normal'
+ NrofInitSamples = 100
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
- MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
- MetaModelOpts.ExpDesign.Method = 'sequential' # 1) normal 2) sequential
+
+ # 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
+ MetaModelOpts.ExpDesign.SamplingMethod = 'hammersley'
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.MaxNSamples = 100
@@ -109,7 +115,7 @@ if __name__ == "__main__":
# -------- Optimality criteria: Optimization ------
# 1)'dual annealing' 2)'BayesOptDesign'
MetaModelOpts.ExpDesign.SeqOptimMethod = 'BayesOptDesign'
- MetaModelOpts.ExpDesign.ExplorationMethod = 'MC' #1)'Voronoi' 2)'LHS' 3) 'MC'
+ MetaModelOpts.ExpDesign.ExplorationMethod = 'random' #1)'Voronoi' 2)'LHS' 3) 'random'
MetaModelOpts.ExpDesign.MaxFunItr = 200
@@ -131,7 +137,7 @@ if __name__ == "__main__":
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Adaptive sparse arbitrary polynomial chaos expansion
- PCEModel = MetaModelOpts.create_PCE(Model)
+ PCEModel = MetaModelOpts.create_metamodel(Model)
#=====================================================
@@ -139,19 +145,16 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 100
- PostPCE.plotFlag = True
-
- PostPCE.ValidMetamodel()
-
- PostPCE.eval_PCEmodel_3D()
+ # Plot to check validation visually.
+ # PostPCE.validMetamodel(nValidSamples=100)
+ # PostPCE.eval_PCEmodel_3D()
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential':
- PostPCE.seqDesignDiagnosticPlots(NrofInitSamples, SaveFig=True)
+ PostPCE.seqDesignDiagnosticPlots()
# Plot the sobol indices
- sobol_cell, total_sobol = PostPCE.SobolIndicesPCE('Z')
+ sobol_cell, total_sobol = PostPCE.sobolIndicesPCE(PlotType='bar')
#=====================================================
#======== Bayesian inference with Emulator ==========
@@ -166,10 +169,11 @@ if __name__ == "__main__":
BayesOpts.PlotPostPred = True
# ----- Define the discrepancy model -------
+ obsData = pd.DataFrame(np.array([5e-1]), columns=Model.Output.Names)
# (Option B)
-# DiscrepancyOpts = Discrepancy('')
-# DiscrepancyOpts.Type = 'Gaussian'
-# DiscrepancyOpts.Parameters = 5e-1
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.Type = 'Gaussian'
+ DiscrepancyOpts.Parameters = obsData**2
# (Option C)
# DiscInputOpts = Input()
@@ -181,7 +185,7 @@ if __name__ == "__main__":
# DiscrepancyOpts = Discrepancy(DiscInputOpts)
-# BayesOpts.Discrepancy = DiscrepancyOpts
+ BayesOpts.Discrepancy = DiscrepancyOpts
Bayes_PCE = BayesOpts.create_Inference()
diff --git a/BayesValidRox/tests/CO2Benchmark/Test_CO2Benchmark.py b/BayesValidRox/tests/CO2Benchmark/Test_CO2Benchmark.py
index 045d4ba1a7686ecf8cba0b9399f1644b71e6cbe6..2608a7a691b3dd9d2083e284203913bb6c257adc 100644
--- a/BayesValidRox/tests/CO2Benchmark/Test_CO2Benchmark.py
+++ b/BayesValidRox/tests/CO2Benchmark/Test_CO2Benchmark.py
@@ -207,24 +207,19 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 3
- PostPCE.plotFlag = True
-
- PostPCE.ValidMetamodel(x_values=np.linspace(0, 8640000, 10)/86400, x_axis="Time [day]")
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=3, x_axis="Time [day]")
- # Comparison with Monte-Carlo reference
- # Compute the moments
- PostPCE.PCEMoments(PCEModel)
+ # Compute the moments and compare with the Monte-Carlo reference
+ PostPCE.plotMoments()
# Plot Mean & Std for all Outputs
- PostPCE.MomentLineplot(PCE_Means = PostPCE.PCEMeans,
- PCE_Stds = PostPCE.PCEStd,
- x_label = 'Time [day]', SaveFig = True)
+ PostPCE.MomentLineplot(x_label = 'Time [day]')
# Plot the evolution of the KLD,BME, and Modified LOOCV error
+ #Information Entropy: -59.51615853207922
if MetaModelOpts.ExpDesign.Method == 'sequential':
- PostPCE.seqDesignDiagnosticPlots(initNSamples, refBME_KLD=(-22.41279, 4.497334360207763), SaveFig=True) #Information Entropy: -59.51615853207922
-
+ PostPCE.seqDesignDiagnosticPlots(refBME_KLD=(-22.41279, 4.497334360207763))
#=====================================================
#============== Bayesian inference =================
diff --git a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
index 1eb4684e8e8dfd57b07e393a9342076c308975f5..1dc65105d5877ce3a8198332fd22c83490d990f7 100755
--- a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
+++ b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
@@ -208,25 +208,16 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 3
- PostPCE.plotFlag = True
- t = np.arange(0, 8640001, 864000)
- PostPCE.ValidMetamodel(x_values=t, x_axis="Time [s]")
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=3)
- # Comparison with Monte-Carlo reference
- # Compute the moments
- PostPCE.PCEMoments(PCEModel)
+ # Compute the moments and compare with the Monte-Carlo reference
+ PostPCE.plotMoments()
- # Plot Mean & Std for all Outputs
- PostPCE.MomentLineplot(PCE_Means = PostPCE.PCEMeans,
- PCE_Stds = PostPCE.PCEStd,
- x_label = 'Time [s]', SaveFig = True)
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential': #infEntropy: -1.1472792737715132
- PostPCE.seqDesignDiagnosticPlots(initNSamples, refBME_KLD=(30.49151, 1.4858325),
- SaveFig=True)
-
+ PostPCE.seqDesignDiagnosticPlots(refBME_KLD=(30.49151, 1.4858325))
#=====================================================
#============== Bayesian inference =================
@@ -320,7 +311,7 @@ if __name__ == "__main__":
ax.plot(value2[xi], value2[yi], "sr")
figPosterior.set_size_inches((10,10))
- figPosterior.savefig('./'+newpath+'/'+key+'.svg',bbox_inches='tight')
+ figPosterior.savefig('./'+newpath+'/'+key+'.pdf',bbox_inches='tight')
theta_true = MetaModelOpts.ExpDesign.MAP
parNames = ['Inj. rate', 'Rel. permeability degree', 'Res. porosity']
diff --git a/BayesValidRox/tests/Fetzer/Test_Fetzer.py b/BayesValidRox/tests/Fetzer/Test_Fetzer.py
index 98454c174038310851ae45b9406ed91eadf22dfc..a7c8c87fa02d7d704787bd8bbd0b6c649aa62d33 100644
--- a/BayesValidRox/tests/Fetzer/Test_Fetzer.py
+++ b/BayesValidRox/tests/Fetzer/Test_Fetzer.py
@@ -22,7 +22,7 @@ except ModuleNotFoundError:
from PyLink.PyLinkForwardModel import PyLinkForwardModel
from surrogate_models.Input import Input
-from surrogate_models.surrogate_models import aPCE
+from surrogate_models.surrogate_models import Metamodel
from PostProcessing.PostProcessing import PostProcessing
from BayesInference.BayesInference import BayesInference, Discrepancy
@@ -86,7 +86,7 @@ if __name__ == "__main__":
#=====================================================
#====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
#=====================================================
- MetaModelOpts = aPCE(Inputs)
+ MetaModelOpts = Metamodel(Inputs)
# Specify the max degree to be compared by the adaptive algorithm:
# The degree with the lowest Leave-One-Out cross-validation (LOO)
@@ -169,7 +169,7 @@ if __name__ == "__main__":
#PCEModel = MetaModelOpts.create_PCE(Model)
MetaModelOpts.slicingforValidation = True
MetaModelOpts.index = 96
- PCEModelCalib, PCEModelValid = MetaModelOpts.create_PCE(Model)
+ PCEModelCalib, PCEModelValid = MetaModelOpts.create_metamodel(Model)
with open('Davarzani2014a_PCEModels.pkl', 'wb') as output:
@@ -181,22 +181,11 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModelCalib)
- ## Compute the moments
- #PostPCE.PCEMoments(PCEModel)
-
- ## Plot Mean & Std for all Outputs
- ## Call the function to create plot
- #PostPCE.MomentLineplot(x_data = PCEModel.ModelOutputDict['Time Steps']
- #, dict_Means = PostPCE.PCEMeans
- #, dict_Stds = PostPCE.PCEStd
- #, x_label = 'Time [day]'
- #, SaveFig = True)
-
- ## Evaluation of surrogate model predictions
- PostPCE.NrofSamples = 1
- PostPCE.plotFlag = True
- t = np.linspace(3600, 691200, 96)/86400
- PostPCE.ValidMetamodel(x_values=t, x_axis="Time [s]")
+ # Compute the moments and compare with the Monte-Carlo reference
+ PostPCE.plotMoments()
+
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=1)
#=====================================================
#======= Bayesian calibration with Emulator =========
diff --git a/BayesValidRox/tests/IshigamiTest/Ishigami.py b/BayesValidRox/tests/IshigamiTest/Ishigami.py
index cd364116a301e96e494e84c4125e7e9b21a3f3e7..a2efe98eed2fe5269ee2e2ce7af59fdfa3b50894 100644
--- a/BayesValidRox/tests/IshigamiTest/Ishigami.py
+++ b/BayesValidRox/tests/IshigamiTest/Ishigami.py
@@ -62,8 +62,7 @@ def Ishigami(xx, *args):
y = term1 + term2 + term3
- return np.vstack((np.zeros(1),y))#np.array(y)[:,None]
- #return np.array([y[:,None]])
+ return np.vstack((np.zeros(1),y))
diff --git a/BayesValidRox/tests/IshigamiTest/Test_Ishigami.py b/BayesValidRox/tests/IshigamiTest/Test_Ishigami.py
index 22bf7d5f595e70225582e1d2596f93cc3a3e1163..c3db37494b58b617353cf4b44930dd1f371b0580 100644
--- a/BayesValidRox/tests/IshigamiTest/Test_Ishigami.py
+++ b/BayesValidRox/tests/IshigamiTest/Test_Ishigami.py
@@ -9,7 +9,7 @@ Created on Wed Jul 10 14:27:35 2019
import numpy as np
-from surrogate_models.surrogate_models import aPCE
+from surrogate_models.surrogate_models import Metamodel
from PyLink.FuncForwardModel import FuncForwardModel
from surrogate_models.Input import Input
from PostProcessing.PostProcessing import PostProcessing
@@ -23,7 +23,7 @@ if __name__ == "__main__":
# TODO: Check if the function is working
# Define model options
- myModel = FuncForwardModel('Ishigami')
+ myModel = FuncForwardModel()
myModel.pyFile = 'Ishigami'
myModel.Type = 'Function'
myModel.pyFile = 'Ishigami'
@@ -55,7 +55,7 @@ if __name__ == "__main__":
#=====================================================
#====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
#=====================================================
- MetaModelOpts = aPCE(Inputs)
+ MetaModelOpts = Metamodel(Inputs)
# Specify the max degree to be compared by the adaptive algorithm:
# The degree with the lowest Leave-One-Out cross-validation (LOO)
@@ -71,8 +71,8 @@ if __name__ == "__main__":
# 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)SGDRegressor: Stochastic gradient descent learning
- MetaModelOpts.RegMethod = 'ARD'
+ # 5)FastARD: Fast Bayesian ARD Regression
+ MetaModelOpts.RegMethod = 'FastARD'
# Print summary of the regression results
#MetaModelOpts.DisplayFlag = True
@@ -83,8 +83,8 @@ if __name__ == "__main__":
MetaModelOpts.addExpDesign()
# One-shot (normal) or Sequential Adaptive (sequential) Design
- MetaModelOpts.ExpDesign.Method = 'sequential'
- NrofInitSamples = 4
+ MetaModelOpts.ExpDesign.Method = 'normal'
+ NrofInitSamples = 100
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
# Sampling methods
@@ -136,27 +136,26 @@ if __name__ == "__main__":
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Adaptive sparse arbitrary polynomial chaos expansion
- PCEModel = MetaModelOpts.create_PCE(myModel)
+ PCEModel = MetaModelOpts.create_metamodel(myModel)
#=====================================================
#========= POST PROCESSING OF METAMODELS ===========
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 1e4
- PostPCE.plotFlag = True
- PostPCE.ValidMetamodel()
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=1e4)
#PostPCE.eval_PCEmodel_3D()
# Compute and print RMSE error
- PostPCE.relErrorPCEModel(nSamples=3000)
+ PostPCE.accuracyCheckMetaModel(nSamples=3000)
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential':
- PostPCE.seqDesignDiagnosticPlots(NrofInitSamples,SaveFig=True)
+ PostPCE.seqDesignDiagnosticPlots()
# Plot the sobol indices
- sobol_cell, total_sobol = PostPCE.SobolIndicesPCE('Z')
+ sobol_cell, total_sobol = PostPCE.sobolIndicesPCE(PlotType='bar')
\ No newline at end of file
diff --git a/BayesValidRox/tests/PollutionTest/Pollution_Test.py b/BayesValidRox/tests/PollutionTest/Pollution_Test.py
index 55d84c1339ff37f1cc299e7d8c49c0614459044d..568885209218e4be97cf2f2a713390f05e67310a 100644
--- a/BayesValidRox/tests/PollutionTest/Pollution_Test.py
+++ b/BayesValidRox/tests/PollutionTest/Pollution_Test.py
@@ -27,7 +27,7 @@ except ModuleNotFoundError:
from PyLink.FuncForwardModel import FuncForwardModel
from surrogate_models.Input import Input
-from surrogate_models.surrogate_models import aPCE
+from surrogate_models.surrogate_models import Metamodel
from PostProcessing.PostProcessing import PostProcessing
from BayesInference.BayesInference import BayesInference, Discrepancy
@@ -45,7 +45,7 @@ if __name__ == "__main__":
Model.Output.Names = ['Z']
# MAP = (8, 0.1, 1, 30.2)
- Model.Observations = np.load("MeasuredData.npy").T
+ Model.Observations['Z'] = np.load("MeasuredData.npy").T
#=====================================================
#========= PROBABILISTIC INPUT MODEL ==============
@@ -81,7 +81,7 @@ if __name__ == "__main__":
#=====================================================
#====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
#=====================================================
- MetaModelOpts = aPCE(Inputs)
+ MetaModelOpts = Metamodel(Inputs)
# Specify the max degree to be compared by the adaptive algorithm:
# The degree with the lowest Leave-One-Out cross-validation (LOO)
@@ -105,10 +105,14 @@ if __name__ == "__main__":
# hypercube sampling of the input model or user-defined values of X and/or Y:
MetaModelOpts.addExpDesign()
- NrofInitSamples = 200
+ # One-shot (normal) or Sequential Adaptive (sequential) Design
+ MetaModelOpts.ExpDesign.Method = 'normal'
+ NrofInitSamples = 100 #75
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
- MetaModelOpts.ExpDesign.SamplingMethod = 'random' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
- MetaModelOpts.ExpDesign.Method = 'normal' # 1) normal 2) sequential
+ # 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
+ MetaModelOpts.ExpDesign.SamplingMethod = 'random'
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.MaxNSamples = 50
@@ -146,7 +150,7 @@ if __name__ == "__main__":
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Adaptive sparse arbitrary polynomial chaos expansion
- PCEModel = MetaModelOpts.create_PCE(Model)
+ PCEModel = MetaModelOpts.create_metamodel(Model)
#=====================================================
@@ -154,15 +158,18 @@ if __name__ == "__main__":
#=====================================================
PostPCE = PostProcessing(PCEModel)
- PostPCE.NrofSamples = 3
- PostPCE.plotFlag = True
-
- #PostPCE.ValidMetamodel()
+ # Plot to check validation visually.
+ PostPCE.validMetamodel(nValidSamples=3)
+
+ # Compute the moments and compare with the Monte-Carlo reference
+ PostPCE.plotMoments()
# Plot the evolution of the KLD,BME, and Modified LOOCV error
if MetaModelOpts.ExpDesign.Method == 'sequential':
- PostPCE.seqDesignDiagnosticPlots(NrofInitSamples, SaveFig=True)
+ PostPCE.seqDesignDiagnosticPlots()
+ import sys
+ sys.exit("STOP")
#=====================================================
#======== Bayesian inference with Emulator ==========
#=====================================================
@@ -195,8 +202,7 @@ if __name__ == "__main__":
Bayes_PCE = BayesOpts.create_Inference()
- import sys
- sys.exit("STOP")
+
# Make directory for saving the plots/outputs
newpath = (r'Outputs_'+MetaModelOpts.ExpDesign.Method+'_'+MetaModelOpts.ExpDesign.SeqOptimMethod)
if not os.path.exists(newpath): os.makedirs(newpath)