From 3ef25b46979c024f9e4ac3b8dc4b8ca012fac4e7 Mon Sep 17 00:00:00 2001
From: faridm69 <faridmohammadi69@gmail.com>
Date: Fri, 24 Jul 2020 17:37:35 +0200
Subject: [PATCH] [all] replaced the eval_PCEmodel functions in all classes
with only the eval_metamodel in metamodel class.
---
.../BayesInference/BayesInference.py | 109 +++------------
BayesValidRox/BayesInference/MCMC.py | 3 +-
.../PostProcessing/PostProcessing.py | 92 +------------
.../surrogate_models/surrogate_models.py | 124 ++++++++++++------
.../Test_AnalyticalFunction.py | 63 ++++-----
5 files changed, 139 insertions(+), 252 deletions(-)
diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 40738d2fb..f48f67a5e 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -88,89 +88,7 @@ class BayesInference:
return self.Samples
- #--------------------------------------------------------------------------------------------------------
- def eval_PCEmodel(self, Samples=None):
-
- PCEModel = self.PCEModel
-
- ModelDict = PCEModel.gp_poly if PCEModel.metaModel.lower() == 'gpe' else PCEModel.CoeffsDict
-
- if Samples is None:
- Samples = self.get_Sample()
- self.Samples = Samples
- else:
- Samples = Samples
- self.NrofSamples = len(Samples)
-
-
- meanPCEOutputs = {}
- stdPCEOutputs = {}
- if PCEModel.metaModel.lower() != 'gpe':
- univ_p_val = PCEModel.univ_basis_vals(Samples)
-
-
- for Outkey, ValuesDict in ModelDict.items():
- PCEOutputs_mean = np.zeros((len(Samples), len(ValuesDict)))
- PCEOutputs_std = np.zeros((len(Samples), len(ValuesDict)))
- idx = 0
- for Inkey, InIdxValues in ValuesDict.items():
-
- # Perdiction with GPE
- if PCEModel.metaModel.lower() == 'gpe':
- gp = PCEModel.gp_poly[Outkey][Inkey]
- y_mean, y_std = gp.predict(Samples, return_std=True)
-
- else: # Perdiction with PCE or pcekriging
- PolynomialDegrees = PCEModel.BasisDict[Outkey][Inkey]
- PSI_Val = PCEModel.PCE_create_Psi(PolynomialDegrees, univ_p_val)
-
- # Perdiction
- if PCEModel.RegMethod.lower() == 'ols': #without error bar
- coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
- y_mean = np.dot(PSI_Val, coeffs)
- y_std = np.zeros(shape=y_mean.shape)
-
- else: #with error bar
- clf_poly = PCEModel.clf_poly[Outkey][Inkey]
- y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
-
- if PCEModel.metaModel == 'PCEKriging':
- gp = PCEModel.gp_poly[Outkey][Inkey]
- y_gp_mean, y_gp_std = gp.predict(Samples, return_std=True)
- y_gp_mean += y_gp_mean
-
- PCEOutputs_mean[:, idx] = y_mean
- PCEOutputs_std[:, idx] = y_std
- idx += 1
-
- if PCEModel.index is None:
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
- stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))
- else:
- meanPCEOutputs[Outkey] = PCEOutputs_mean
- stdPCEOutputs[Outkey] = PCEOutputs_std
- else:
- index = PCEModel.index
- if self.Name == 'Calib':
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- meanPCEOutputs[Outkey] = PCA.mean_[:index] + np.dot(PCEOutputs_mean,PCA.components_)[:,:index]
- stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,:index]
- else:
- meanPCEOutputs[Outkey] = PCEOutputs_mean[:,:index]
- stdPCEOutputs[Outkey] = PCEOutputs_std[:,:index]
- else:
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- meanPCEOutputs[Outkey] = PCA.mean_[index:] + np.dot(PCEOutputs_mean,PCA.components_)[:,index:]
- stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,index:]
- else:
- meanPCEOutputs[Outkey] = PCEOutputs_mean[:,index:]
- stdPCEOutputs[Outkey] = PCEOutputs_std[:,index:]
-
- return meanPCEOutputs, stdPCEOutputs
+
#--------------------------------------------------------------------------------------------------------
def eval_Model(self, Samples=None):
@@ -319,7 +237,8 @@ class BayesInference:
"""
Calculates the correction factor for BMEs.
"""
- OrigModelOutput = self.PCEModel.ExpDesign.Y
+ PCEModel = self.PCEModel
+ OrigModelOutput = PCEModel.ExpDesign.Y
# Extract the requested model outputs for likelihood calulation
@@ -338,8 +257,8 @@ class BayesInference:
NrofBayesSamples = self.NrofSamples
# Evaluate PCEModel on the experimental design
- Samples = self.PCEModel.ExpDesign.X
- OutputRS, stdOutputRS = self.eval_PCEmodel(Samples)
+ Samples = PCEModel.ExpDesign.X
+ OutputRS, stdOutputRS = PCEModel.eval_metamodel(samples=Samples)
# Reset the NrofSamples to NrofBayesSamples
self.NrofSamples = NrofBayesSamples
@@ -463,7 +382,8 @@ class BayesInference:
#--------------------------------------------------------------------------------------------------------
def PosteriorPredictive(self):
- Model = self.PCEModel.ModelObj
+ PCEModel = self.PCEModel
+ Model = PCEModel.ModelObj
# Make a directory to save the prior/posterior predictive
OutputDir = (r'Outputs_Bayes_' + Model.Name + '_' + self.Name)
@@ -495,7 +415,7 @@ class BayesInference:
# Prior predictive
if self.emulator:
PriorPred = self.meanPCEPriorPred
- PosteriorPred, _ = self.eval_PCEmodel(Samples=Posterior_df.to_numpy())
+ PosteriorPred, _ = PCEModel.eval_metamodel(samples=Posterior_df.to_numpy())
else:
PriorPred = self.ModelPriorPred
PosteriorPred = self.eval_Model(Samples=Posterior_df.to_numpy())
@@ -534,14 +454,17 @@ class BayesInference:
def create_Inference(self):
- Model = self.PCEModel.ModelObj
- NofPa = self.PCEModel.NofPa
+ PCEModel = self.PCEModel
+ Model = PCEModel.ModelObj
+ NofPa = PCEModel.NofPa
OutputNames = Model.Output.Names
BootstrapItrNr = self.BootstrapItrNr
# If the prior is set by the user, take it.
- if self.Samples is not None:
+ if self.Samples is None:
+ self.Samples = self.get_Sample()
+ else:
try:
SamplesAllParameters = self.Samples.to_numpy()
except:
@@ -636,7 +559,7 @@ class BayesInference:
# ---------------- Likelihood calculation ----------------
if self.emulator:
# Evaluate the PCEModel
- self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
+ self.meanPCEPriorPred, self.stdPCEPriorPred = PCEModel.eval_metamodel(samples=self.Samples)
# unknown sigma2
if optSigma == 'C':
@@ -699,7 +622,7 @@ class BayesInference:
if self.PlotPostPred:
# Run the models for MAP
# PCEModel
- MAP_PCEModel, MAP_PCEModelstd = self.eval_PCEmodel(Samples=MAP_theta)
+ MAP_PCEModel, MAP_PCEModelstd = PCEModel.eval_metamodel(samples=MAP_theta)
self.MAPpceModelMean = MAP_PCEModel
self.MAPpceModelStd = MAP_PCEModelstd
diff --git a/BayesValidRox/BayesInference/MCMC.py b/BayesValidRox/BayesInference/MCMC.py
index 42e875bca..db1c4a4a9 100755
--- a/BayesValidRox/BayesInference/MCMC.py
+++ b/BayesValidRox/BayesInference/MCMC.py
@@ -123,12 +123,13 @@ class MCMC():
def log_likelihood(self, theta, Observation, TotalSigma2):
BayesOpts = self.BayesOpts
+ PCEModel = BayesOpts.PCEModel
ndimTheta = theta.ndim
theta = theta if ndimTheta != 1 else theta.reshape((1,-1))
if BayesOpts.emulator:
# Evaluate the PCEModel
- meanPCEPriorPred, BayesOpts.stdPCEPriorPred = BayesOpts.eval_PCEmodel(theta)
+ meanPCEPriorPred, BayesOpts.stdPCEPriorPred = PCEModel.eval_metamodel(samples=theta)
else:
# Evaluate the origModel
meanPCEPriorPred = BayesOpts.eval_Model(theta)
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 9899bde07..3fa0d9ac4 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -132,91 +132,6 @@ class PostProcessing:
return self.Samples
- #--------------------------------------------------------------------------------------------------------
- def eval_PCEmodel(self, Samples=None):
-
- PCEModel = self.PCEModel
-
- ModelDict = PCEModel.gp_poly if PCEModel.metaModel.lower() == 'gpe' else PCEModel.CoeffsDict
-
- if Samples is None:
- Samples = self.get_Sample()
- self.Samples = Samples
- else:
- Samples = Samples
- self.NrofSamples = len(Samples)
-
- self.PCEOutputs = {}
- self.PCEOutputs_std = {}
- if PCEModel.metaModel.lower() != 'gpe':
- univ_p_val = PCEModel.univ_basis_vals(Samples)
-
- # Loop over outputs
- for Outkey, ValuesDict in ModelDict.items():
-
- PCEOutputs_mean = np.zeros((len(Samples), len(ValuesDict)))
- PCEOutputs_std = np.zeros((len(Samples), len(ValuesDict)))
-
- idx = 0
- for Inkey, InIdxValues in ValuesDict.items():
-
- # Perdiction with GPE
- if PCEModel.metaModel.lower() == 'gpe':
- gp = PCEModel.gp_poly[Outkey][Inkey]
- y_mean, y_std = gp.predict(Samples, return_std=True)
-
- else: # Perdiction with PCE or pcekriging
- PolynomialDegrees = PCEModel.BasisDict[Outkey][Inkey]
- clf_poly = PCEModel.clf_poly[Outkey][Inkey]
- PSI_Val = PCEModel.PCE_create_Psi(PolynomialDegrees, univ_p_val)
- try:
- # with error bar
- y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
- except:
- # without error bar
- coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
- y_mean = np.dot(PSI_Val, coeffs)
- y_std = np.zeros(shape=y_mean.shape)
-
- if PCEModel.metaModel.lower() == 'pcekriging':
- gp = PCEModel.gp_poly[Outkey][Inkey]
- y_gp_mean, y_gp_std = gp.predict(Samples, return_std=True)
- y_mean += y_gp_mean
- y_std = y_gp_std
-
- PCEOutputs_mean[:, idx] = y_mean
- PCEOutputs_std[:, idx] = y_std
- idx += 1
-
- if PCEModel.index is None:
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- self.PCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
- self.PCEOutputs_std[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))
- else:
- self.PCEOutputs[Outkey] = PCEOutputs_mean
- self.PCEOutputs_std[Outkey] = PCEOutputs_std
- else:
- index = PCEModel.index
- if self.Name.lower() == 'calib':
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- self.PCEOutputs[Outkey] = PCA.mean_[:index] + np.dot(PCEOutputs_mean,PCA.components_)[:,:index]
- self.PCEOutputs_std[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,:index]
- else:
- self.PCEOutputs[Outkey] = PCEOutputs_mean[:,:index]
- self.PCEOutputs_std[Outkey] = PCEOutputs_std[:,:index]
- else:
- if PCEModel.DimRedMethod.lower() == 'pca':
- PCA = PCEModel.pca[Outkey]
- self.PCEOutputs[Outkey] = PCA.mean_[index:] + np.dot(PCEOutputs_mean,PCA.components_)[:,index:]
- self.PCEOutputs_std[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,index:]
- else:
- self.PCEOutputs[Outkey] = PCEOutputs_mean[:,index:]
- self.PCEOutputs_std[Outkey] = PCEOutputs_std[:,index:]
-
- return self.PCEOutputs, self.PCEOutputs_std
-
#--------------------------------------------------------------------------------------------------------
def eval_PCEmodel_3D(self, SaveFig=True):
@@ -353,10 +268,11 @@ class PostProcessing:
"""
Evaluate the meta model and the PCEModel
"""
+ metaModel = self.PCEModel
Samples = self.get_Sample()
self.eval_Model(Samples)
- self.eval_PCEmodel(Samples)
+ self.PCEOutputs, self.PCEOutputs_std = metaModel.eval_metamodel(samples=Samples)
if self.plotFlag == True:
try:
@@ -377,6 +293,8 @@ class PostProcessing:
"""
Evaluate the relative error of the PCEModel
"""
+ metaModel = self.PCEModel
+
# Set the number of samples
self.NrofSamples = nSamples if nSamples is not None else Samples.shape[0]
@@ -387,7 +305,7 @@ class PostProcessing:
ModelOutputs = self.eval_Model(Samples) if validOutputsDict is None else validOutputsDict
# Run the PCE model with the generated samples
- PCEOutputs, PCEOutputs_std = self.eval_PCEmodel(Samples)
+ PCEOutputs, PCEOutputs_std = metaModel.eval_metamodel(samples=Samples)
self.RMSE = {}
# Loop over the keys and compute RMSE error.
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 05cf52369..4431814f0 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -166,6 +166,7 @@ class Metamodel:
return self.PolynomialDegrees
#--------------------------------------------------------------------------------------------------------
+
def addExpDesign(self):
self.ExpDesign = ExpDesigns(self.Inputs)
@@ -1120,7 +1121,7 @@ class Metamodel:
if UtilMethod == 'Entropy':
# ----- Entropy/MMSE/active learning MacKay(ALM) -----
# Compute perdiction variance of the old model
- Y_PC_can, std_PC_can = self.eval_PCEmodel(X_can)
+ Y_PC_can, std_PC_can = self.eval_metamodel(samples=X_can)
canPredVar = {key:std_PC_can[key]**2 for key in OutputNames}
varPCE = []
@@ -1131,7 +1132,7 @@ class Metamodel:
elif UtilMethod == 'EIGF':
# ----- Expected Improvement for Global fit -----
# Compute perdiction error and variance of the old model
- Y_PC_can, std_PC_can = self.eval_PCEmodel(X_can)
+ Y_PC_can, std_PC_can = self.eval_metamodel(samples=X_can)
predError = {key:Y_PC_can[key] for key in OutputNames}
canPredVar = {key:std_PC_can[key]**2 for key in OutputNames}
@@ -1150,7 +1151,7 @@ class Metamodel:
nTriningPoints = self.ExpDesign.X.shape[0]
# Evaluate the PCE metamodels at that location ???
- Y_mean_can, Y_std_can = self.eval_PCEmodel(np.array([X_can]))
+ Y_mean_can, Y_std_can = self.eval_metamodel(samples=np.array([X_can]))
# Get the data
ObservationData = self.Observations
@@ -1305,7 +1306,7 @@ class Metamodel:
oldExpDesignY = PCEModel.ExpDesign.Y
# Evaluate the PCE metamodels at that location ???
- Y_PC_can, _ = self.eval_PCEmodel(np.array([X_can]))
+ Y_PC_can, _ = self.eval_metamodel(samples=np.array([X_can]))
# Add all suggestion as new ExpDesign
NewExpDesignX = np.vstack((oldExpDesignX, X_can))
@@ -1339,7 +1340,7 @@ class Metamodel:
X_MC = self.ExpDesign.GetSample(MCsize, 'random');
# Evaluate the PCEModel at the given samples
- Y_PC_MC, std_PC_MC = PCE_SparseBayes_can.eval_PCEmodel(X_MC)
+ Y_PC_MC, std_PC_MC = PCE_SparseBayes_can.eval_metamodel(samples=X_MC)
# Likelihood computation (Comparison of data
# and simulation results via PCE with candidate design)
@@ -1541,7 +1542,7 @@ class Metamodel:
# TODO: New adaptive trade-off according to Liu et al. (2017)
# Mean squared error for last design point
- lastPCEY, _ = self.eval_PCEmodel(OldExpDesign[-1].reshape(1,-1))
+ lastPCEY, _ = self.eval_metamodel(samples=OldExpDesign[-1].reshape(1,-1))
if 'x_values' in OutputDictY: del OutputDictY['x_values']
mseError = mean_squared_error(np.array(list(lastPCEY.values()))[:,0], np.array(list(OutputDictY.values()))[:,-1,:])
@@ -2002,52 +2003,95 @@ class Metamodel:
return Xnew
#--------------------------------------------------------------------------------------------------------
- def eval_PCEmodel(self,Samples):
+ def eval_metamodel(self,samples=None, nsamples=None, samplingMethod='random', return_samples=False):
+ ModelDict = self.gp_poly if self.metaModel.lower() == 'gpe' else self.CoeffsDict
- univ_p_val = self.univ_basis_vals(Samples)
+ if samples is None:
+ if nsamples is None:
+ self.NrofSamples = 100000
+ else:
+ self.NrofSamples = nsamples
+
+ self.Samples = self.ExpDesign.GetSample(self.NrofSamples, samplingMethod)
+ else:
+ self.Samples = samples
+ self.NrofSamples = len(samples)
+
+ if self.metaModel.lower() != 'gpe':
+ univ_p_val = self.univ_basis_vals(self.Samples)
meanPCEOutputs = {}
stdPCEOutputs = {}
- for Outkey, ValuesDict in self.CoeffsDict.items():
+ # Loop over outputs
+ for Outkey, ValuesDict in ModelDict.items():
- PCEOutputs_mean = np.zeros((len(Samples), len(ValuesDict)))
- PCEOutputs_std = np.zeros((len(Samples), len(ValuesDict)))
+ PCEOutputs_mean = np.zeros((len(self.Samples), len(ValuesDict)))
+ PCEOutputs_std = np.zeros((len(self.Samples), len(ValuesDict)))
+ idx = 0
for Inkey, InIdxValues in ValuesDict.items():
- idx = int(Inkey.split('_')[1]) - 1
- PolynomialDegrees = self.BasisDict[Outkey][Inkey]
- PSI_Val = self.PCE_create_Psi(PolynomialDegrees, univ_p_val)
-
- # Perdiction
- try: # with error bar
- clf_poly = self.clf_poly[Outkey][Inkey]
- y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
-
- except: # without error bar
- coeffs = self.CoeffsDict[Outkey][Inkey]
- y_mean = np.dot(PSI_Val, coeffs)
- y_std = np.zeros(shape=y_mean.shape)
-
- if self.metaModel.lower() == 'pcekriging':
+ # Perdiction with GPE
+ if self.metaModel.lower() == 'gpe':
gp = self.gp_poly[Outkey][Inkey]
- y_gp_mean, y_gp_std = gp.predict(Samples, return_std=True)
- y_mean += y_gp_mean
- y_std = y_gp_std
+ y_mean, y_std = gp.predict(self.Samples, return_std=True)
+
+ else: # Perdiction with PCE or pcekriging
+ PolynomialDegrees = self.BasisDict[Outkey][Inkey]
+ PSI_Val = self.PCE_create_Psi(PolynomialDegrees, univ_p_val)
+
+ # Perdiction
+ try: # with error bar
+ clf_poly = self.clf_poly[Outkey][Inkey]
+ y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
+
+ except: # without error bar
+ coeffs = self.CoeffsDict[Outkey][Inkey]
+ y_mean = np.dot(PSI_Val, coeffs)
+ y_std = np.zeros(shape=y_mean.shape)
+
+ if self.metaModel.lower() == 'pcekriging':
+ gp = self.gp_poly[Outkey][Inkey]
+ y_gp_mean, y_gp_std = gp.predict(self.Samples, return_std=True)
+ y_mean += y_gp_mean
+ y_std = y_gp_std
PCEOutputs_mean[:, idx] = y_mean
PCEOutputs_std[:, idx] = y_std
+ idx += 1
- if self.DimRedMethod.lower() == 'pca':
- PCA = self.pca[Outkey]
- meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
- stdPCEOutputs[Outkey] = np.dot(PCEOutputs_std, PCA.components_)
+ if self.index is None:
+ if self.DimRedMethod.lower() == 'pca':
+ PCA = self.pca[Outkey]
+ meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_)
+ stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))
+ else:
+ meanPCEOutputs[Outkey] = PCEOutputs_mean
+ stdPCEOutputs[Outkey] = PCEOutputs_std
else:
- meanPCEOutputs[Outkey] = PCEOutputs_mean
- stdPCEOutputs[Outkey] = PCEOutputs_std
-
- return meanPCEOutputs, stdPCEOutputs
+ index = self.index
+ if self.Name.lower() == 'calib':
+ if self.DimRedMethod.lower() == 'pca':
+ PCA = self.pca[Outkey]
+ meanPCEOutputs[Outkey] = PCA.mean_[:index] + np.dot(PCEOutputs_mean,PCA.components_)[:,:index]
+ stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,:index]
+ else:
+ meanPCEOutputs[Outkey] = PCEOutputs_mean[:,:index]
+ stdPCEOutputs[Outkey] = PCEOutputs_std[:,:index]
+ else:
+ if self.DimRedMethod.lower() == 'pca':
+ PCA = self.pca[Outkey]
+ meanPCEOutputs[Outkey] = PCA.mean_[index:] + np.dot(PCEOutputs_mean,PCA.components_)[:,index:]
+ stdPCEOutputs[Outkey] = np.sqrt(np.dot(PCEOutputs_std**2, PCA.components_**2))[:,index:]
+ else:
+ meanPCEOutputs[Outkey] = PCEOutputs_mean[:,index:]
+ stdPCEOutputs[Outkey] = PCEOutputs_std[:,index:]
+
+ if return_samples:
+ return meanPCEOutputs, stdPCEOutputs, self.Samples
+ else:
+ return meanPCEOutputs, stdPCEOutputs
#--------------------------------------------------------------------------------------------------------
def normpdf(self, PCEOutputs, std_PC_MC, ObservationData, Sigma2s):
@@ -2109,7 +2153,7 @@ class Metamodel:
# ----- PCE model outputs on the ExpDesign ---------------
# Evaluate PCEModel on the experimental design
Samples = PCEModel.ExpDesign.X
- OutputRS, stdOutputRS = PCEModel.eval_PCEmodel(Samples)
+ OutputRS, stdOutputRS = PCEModel.eval_metamodel(samples=Samples)
# Flatten the mean Outputs for PCEModel
TotalPCEOutputs = np.empty((SampleSize,0))
@@ -2324,7 +2368,7 @@ class Metamodel:
MCsize = X_MC.shape[0]
# Monte Carlo simulation for the candidate design
- Y_PC_MC, std_PC_MC = PCEModel.eval_PCEmodel(X_MC)
+ Y_PC_MC, std_PC_MC = PCEModel.eval_metamodel(samples=X_MC)
# Likelihood computation (Comparison of data and
# simulation results via PCE with candidate design)
@@ -2410,7 +2454,7 @@ class Metamodel:
ModelOutputs = self.validModelRuns
# Run the PCE model with the generated samples
- PCEOutputs, PCEOutputs_std = self.eval_PCEmodel(Samples)
+ PCEOutputs, PCEOutputs_std = self.eval_metamodel(samples=Samples)
validError_dict = {}
# Loop over the keys and compute RMSE error.
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index 26a2c0ade..af95523e8 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -70,18 +70,18 @@ if __name__ == "__main__":
# standard deviation
Inputs = Input()
- 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]
+ # 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 ============
@@ -127,13 +127,13 @@ if __name__ == "__main__":
# One-shot (normal) or Sequential Adaptive (sequential) Design
MetaModelOpts.ExpDesign.Method = 'normal'
- NrofInitSamples = 100 #75
+ NrofInitSamples = 150 #75
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
- MetaModelOpts.ExpDesign.SamplingMethod = 'latin_hypercube'
+ MetaModelOpts.ExpDesign.SamplingMethod = 'random'
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
@@ -208,29 +208,30 @@ if __name__ == "__main__":
t = np.arange(0, 10, 1.) / 9
PostPCE.ValidMetamodel(x_values=t, x_axis="Time [s]")
- # # Comparison with Monte-Carlo reference
- # # Compute the moments
- # PostPCE.PCEMoments(PCEModel)
+ # 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)
+ # 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':
- # refBME_KLD = np.load("data/refBME_KLD_"+str(ndim)+".npy")
- # PostPCE.seqDesignDiagnosticPlots(NrofInitSamples, refBME_KLD
- # , SaveFig=True)
+ # 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)
- # # Plot the sobol indices
- # sobol_cell, total_sobol = PostPCE.SobolIndicesPCE(x_values=t, x_axis="Time [s]",
- # SaveFig=True)
+ # Plot the sobol indices
+ sobol_cell, total_sobol = PostPCE.SobolIndicesPCE(x_values=t, x_axis="Time [s]",
+ SaveFig=True)
# Compute and print RMSE error
PostPCE.relErrorPCEModel(nSamples=3000)
- sys.exit('STOP!!')
+
+ # sys.exit('STOP!!')
#=====================================================
#======== Bayesian inference with Emulator ==========
#=====================================================
@@ -238,14 +239,14 @@ if __name__ == "__main__":
BayesOpts.emulator = True
# BME Bootstrap
- BayesOpts.Bootstrap = True
+ BayesOpts.Bootstrap = False
BayesOpts.BootstrapItrNr = 10
BayesOpts.BootstrapNoise = 0.05
# Select the inference method
- # BayesOpts.SamplingMethod = 'MCMC'
- # BayesOpts.MCMCnSteps = 1000
- # BayesOpts.MultiProcessMCMC = True
+ BayesOpts.SamplingMethod = 'MCMC'
+ BayesOpts.MCMCnSteps = 1000
+ BayesOpts.MultiProcessMCMC = True
BayesOpts.PlotPostPred = True
--
GitLab