From bc36444f4e4bc1cf613df67c01afd9c34530e28b Mon Sep 17 00:00:00 2001
From: faridm69 <faridmohammadi69@gmail.com>
Date: Wed, 17 Jun 2020 13:32:59 +0200
Subject: [PATCH] [surrogate] added principal component analysis.
---
.../BayesInference/BayesInference.py | 453 +++++++++++-------
.../PostProcessing/PostProcessing.py | 59 +--
.../surrogate_models/RegressionFastARD.py | 2 +-
.../surrogate_models/surrogate_models.py | 236 +++++----
.../AnalytFuncValid_Test.py | 11 +-
.../Test_AnalyticalFunction.py | 18 +-
6 files changed, 459 insertions(+), 320 deletions(-)
diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index e734a1319..5c6032abc 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -27,12 +27,14 @@ plt.rc('savefig', dpi=1000)
from .Discrepancy import Discrepancy
+from .MCMC import MCMC
from surrogate_models.ExpDesigns import ExpDesigns
class BayesInference:
def __init__(self, PCEModel):
self.PCEModel = PCEModel
self.Name = 'Calib'
+ self.SamplingMethod = 'rejection'
self.MAP ='Mean'
self.PCEMeans = {}
self.PCEStd = {}
@@ -80,7 +82,7 @@ class BayesInference:
NrSamples = self.NrofSamples
- self.Samples = PCEModel.ExpDesign.GetSample(NrSamples, SamplingMethod='random')
+ self.Samples = PCEModel.ExpDesign.GetSample(NrSamples, 'random')
return self.Samples
@@ -120,11 +122,28 @@ class BayesInference:
coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
PCEOutputs_mean[:, idx] = np.dot(PSI_Val, coeffs)
+ if PCEModel.metaModel == '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
+
idx += 1
-
- meanPCEOutputs[Outkey] = PCEOutputs_mean
- stdPCEOutputs[Outkey] = PCEOutputs_std
+ PCA = PCEModel.pca[Outkey]
+ if PCEModel.index is None:
+ meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_) # PCEOutputs_mean
+ stdPCEOutputs[Outkey] = np.dot(PCEOutputs_std, PCA.components_) # PCEOutputs_std
+ else:
+ index = PCEModel.index
+ if self.Name == 'Calib':
+ meanPCEOutputs[Outkey] = PCA.mean_[:index] + np.dot(PCEOutputs_mean,PCA.components_)[:,:index] # PCEOutputs_mean
+ stdPCEOutputs[Outkey] = np.dot(PCEOutputs_std, PCA.components_)[:,:index] # PCEOutputs_std
+ else:
+ meanPCEOutputs[Outkey] = PCA.mean_[index:] + np.dot(PCEOutputs_mean,PCA.components_)[:,index:] # PCEOutputs_mean
+ stdPCEOutputs[Outkey] = np.dot(PCEOutputs_std, PCA.components_)[:,index:] # PCEOutputs_std
+
+
return meanPCEOutputs, stdPCEOutputs
#--------------------------------------------------------------------------------------------------------
@@ -197,7 +216,7 @@ class BayesInference:
meanPCE = np.empty((SampleSize,0))
stdPCE = np.empty((SampleSize,0))
for outputType in OutputType:
- meanPCE = np.hstack((meanPCE, self.meanPCEPriorPred[outputType]))
+ meanPCE = np.hstack((meanPCE, Outputs[outputType]))#self.meanPCEPriorPred[outputType]))
stdPCE = np.hstack((stdPCE, self.stdPCEPriorPred[outputType]))
# Expected value of variance (Assump: i.i.d stds)
@@ -419,7 +438,6 @@ class BayesInference:
return self.Posterior_df
-
#--------------------------------------------------------------------------------------------------------
def PosteriorPredictive(self):
@@ -434,7 +452,10 @@ class BayesInference:
MeasuredData = Model.Observations
else:
MeasuredData = Model.read_ObservationValid()
-
+
+ # X_values
+ x_key = list(MeasuredData)[0]
+
try:
Sigma2Prior = self.Discrepancy.Sigma2Prior
except:
@@ -446,35 +467,31 @@ class BayesInference:
if Sigma2Prior is not None:
Posterior_df = Posterior_df.drop(labels=self.Discrepancy.Name, axis=1)
else:
- Posterior_df =self.Posterior_df
-
-
+ Posterior_df = self.Posterior_df
- if self.emulator == True:
+ # Prior predictive
+ if self.emulator:
PriorPred = self.meanPCEPriorPred
PosteriorPred, _ = self.eval_PCEmodel(Samples=Posterior_df.to_numpy())
-
else:
-
PriorPred = self.ModelPriorPred
PosteriorPred = self.eval_Model(Samples=Posterior_df.to_numpy())
- # Convert the array to a categorical data frame (PriorPredictive)
- # X_values
- x_key = list(MeasuredData)[0]
-
- X_values = np.repeat(MeasuredData[x_key] , PriorPred[Model.Output.Names[0]].shape[0])
-
+ self.PCEAllPred_df = PosteriorPred
- # Prior Predictive
- priorh5File = "priorPredictive.h5"
- for Out_idx, OutputName in enumerate(Model.Output.Names):
- PriorPred_dict = {}
- PriorPred_dict[x_key] = X_values
- PriorPred_dict[OutputName] = PriorPred[OutputName].flatten('F')
- PriorPred_df = pd.DataFrame(PriorPred_dict)
- PriorPred_df.to_hdf(OutputDir+'/'+priorh5File, "/data/"+OutputName)
+ if self.SamplingMethod == 'rejection':
+ # Convert the array to a categorical data frame (PriorPredictive)
+ X_values = np.repeat(MeasuredData[x_key] , PriorPred[Model.Output.Names[0]].shape[0])
+
+ # create the dataframe
+ priorh5File = "priorPredictive.h5"
+ for Out_idx, OutputName in enumerate(Model.Output.Names):
+ PriorPred_dict = {}
+ PriorPred_dict[x_key] = X_values
+ PriorPred_dict[OutputName] = PriorPred[OutputName].flatten('F')
+ PriorPred_df = pd.DataFrame(PriorPred_dict)
+ PriorPred_df.to_hdf(OutputDir+'/'+priorh5File, "/data/"+OutputName)
# Convert the array to a categorical data frame (PosteriorPredictive)
# X_values
@@ -575,7 +592,7 @@ class BayesInference:
# ---------------- Bootstrap & TOM --------------------
- if self.Bootstrap is True:
+ if self.Bootstrap:
if len(self.perturbedData) != 0:
Data = self.perturbedData
else:
@@ -593,93 +610,64 @@ class BayesInference:
# -----------------------------------------------------
- # Loop over the perturbed observation data
+ # ----- Loop over the perturbed observation data ------
# -----------------------------------------------------
- for itrIdx, data in tqdm(enumerate(Data), ascii=True, desc ="Boostraping the BME calculations"):
-
- # ---------------- Likelihood calculation ----------------
- if self.emulator == True:
- # Evaluate the PCEModel
- self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
+ if self.SamplingMethod == 'rejection':
+ for itrIdx, data in tqdm(enumerate(Data), ascii=True, desc ="Boostraping the BME calculations"):
- # unknown sigma2
- if optSigma == 'C':
- Likelihoods[:,itrIdx] = self.normpdfSigma2(self.meanPCEPriorPred, data, TotalSigma2)
- else:
- # known sigma2
- Likelihoods[:,itrIdx] = self.normpdf(self.meanPCEPriorPred, data, TotalSigma2)
-
- else:
- self.ModelPriorPred = self.eval_Model(Samples=self.Samples)
- # unknown sigma2
- if optSigma == 'C':
- Likelihoods[:,itrIdx] = self.normpdfSigma2(self.ModelPriorPred, data, TotalSigma2)
- else:
- # known sigma2
- Likelihoods[:,itrIdx] = self.normpdf(self.ModelPriorPred, data, TotalSigma2)
-
- # ---------------- BME Calculations ----------------
- # BME (log)
- BME[itrIdx] = np.log(np.nanmean(np.exp(Likelihoods[:,itrIdx])))
-
- # BME correction when using Emulator
- # if self.emulator:
- # BME_Corr[itrIdx] = self.BME_Corr_Weight(data, sigma2)
-
- # ---------------- Store BME, Likelihoods for all ----------------
- # Likelihoods (Size: NrofSamples,BootstrapItrNr)
- self.Likelihoods = Likelihoods
-
- # BME (log) (Size: 1,BootstrapItrNr)
- self.BME = BME
-
- # TODO: BMECorrFactor (log) (Size: 1,BootstrapItrNr)
- #if self.emulator: self.BMECorrFactor = BME_Corr
-
- # BME = BME + BMECorrFactor
- if self.emulator: self.BME = self.BME #+ self.BMECorrFactor
-
- # ---------------- Rejection Sampling ----------------
- self.Posterior_df = self.Rejection_Sampling()
-
- # Posterior perdictive for PCEModel
- AllPred_df = self.PosteriorPredictive()
-
- # ---------------- Plots ----------------
- # Plot posteior parameters
- OutputDir = (r'Outputs_Bayes_' + Model.Name + '_' + self.Name)
- if not os.path.exists(OutputDir): os.makedirs(OutputDir)
+ # ---------------- Likelihood calculation ----------------
+ if self.emulator:
+ # Evaluate the PCEModel
+ self.meanPCEPriorPred, self.stdPCEPriorPred = self.eval_PCEmodel(Samples=self.Samples)
+
+ # unknown sigma2
+ if optSigma == 'C':
+ Likelihoods[:,itrIdx] = self.normpdfSigma2(self.meanPCEPriorPred, data, TotalSigma2)
+ else:
+ # known sigma2
+ Likelihoods[:,itrIdx] = self.normpdf(self.meanPCEPriorPred, data, TotalSigma2)
- # Plot logBME
- if self.BootstrapItrNr != 1:
- # Computing the TOM performance
- x_values = np.linspace(np.min(self.BME), np.max(self.BME), 1000)
- dof = ObservationData.shape[0]
- self.logTOMBME = stats.chi2(dof).logpdf(-1*x_values)#pdf(x_values)
-
- fig, ax = plt.subplots()
- sns.kdeplot(self.logTOMBME, ax=ax, color="green", shade=True)
- sns.kdeplot(self.BME, ax=ax, color="blue", shade=True,label='Model BME')
-
- ax.set_xlabel('log$_{10}$(BME)')
- ax.set_ylabel('Probability density')
+ else:
+ self.ModelPriorPred = self.eval_Model(Samples=self.Samples)
+ # unknown sigma2
+ if optSigma == 'C':
+ Likelihoods[:,itrIdx] = self.normpdfSigma2(self.ModelPriorPred, data, TotalSigma2)
+ else:
+ # known sigma2
+ Likelihoods[:,itrIdx] = self.normpdf(self.ModelPriorPred, data, TotalSigma2)
+
+ # ---------------- BME Calculations ----------------
+ # BME (log)
+ BME[itrIdx] = np.log(np.nanmean(np.exp(Likelihoods[:,itrIdx])))
+
+ # BME correction when using Emulator
+ # if self.emulator:
+ # BME_Corr[itrIdx] = self.BME_Corr_Weight(data, sigma2)
+
+ # ---------------- Store BME, Likelihoods for all ----------------
+ # Likelihoods (Size: NrofSamples,BootstrapItrNr)
+ self.Likelihoods = Likelihoods
- from matplotlib.patches import Patch
- legend_elements = [Patch(facecolor='green', edgecolor='green',label='TOM BME'),
- Patch(facecolor='blue', edgecolor='blue',label='Model BME')]
- ax.legend(handles=legend_elements,fontsize=14)
+ # BME (log) (Size: 1,BootstrapItrNr)
+ self.BME = BME
+
+ # TODO: BMECorrFactor (log) (Size: 1,BootstrapItrNr)
+ #if self.emulator: self.BMECorrFactor = BME_Corr
+
+ # BME = BME + BMECorrFactor
+ if self.emulator: self.BME = self.BME #+ self.BMECorrFactor
- if self.emulator:
- plotname = '/BME_hist_'+Model.Name+'_emulator'
- else:
- plotname = '/BME_hist_'+Model.Name
+ # ---------------- Rejection Sampling ----------------
+ self.Posterior_df = self.Rejection_Sampling()
+ else:
+ # TODO: MCMC
+ initsamples = None if self.Samples is None else self.Samples
+ nsteps = 1000 #if self.NrofSamples is None else self.NrofSamples
+ MCMC_ = MCMC(self, initsamples=initsamples, nwalkers=2*self.PCEModel.NofPa,
+ nsteps = nsteps)
+ self.Posterior_df = MCMC_.run_sampler(Data[0], TotalSigma2)
- plt.savefig('./'+OutputDir+plotname+'.svg', bbox_inches = 'tight')
- plt.show()
- plt.close()
-
-
# -------- Find MAP and run PCEModel and origModel --------
# Compute the MAP
if self.MAP =='Mean':
@@ -687,19 +675,31 @@ class BayesInference:
else:
MAP_theta = stats.mode(self.Posterior_df.to_numpy(),axis=0)[0]
print("\nPoint estimator:\n", MAP_theta[0])
-
- # PCEModel
- MAP_PCEModel, MAP_PCEModelstd = self.eval_PCEmodel(Samples=MAP_theta)
- self.MAPpceModelMean = MAP_PCEModel
- self.MAPpceModelStd = MAP_PCEModelstd
- # origModel
- MAP_origModel = self.eval_Model(Samples=MAP_theta)
- self.MAPorigModel = MAP_origModel
+ if self.PlotPostPred:
+ # Run the models for MAP
+ # PCEModel
+ MAP_PCEModel, MAP_PCEModelstd = self.eval_PCEmodel(Samples=MAP_theta)
+ self.MAPpceModelMean = MAP_PCEModel
+ self.MAPpceModelStd = MAP_PCEModelstd
+
+ # origModel
+ MAP_origModel = self.eval_Model(Samples=MAP_theta)
+ self.MAPorigModel = MAP_origModel
+
+ # -------- Posterior perdictive -----------
+ self.PosteriorPredictive()
+
+ # -----------------------------------------------------
+ # ------------------ Visualization --------------------
+ # -----------------------------------------------------
# -------- Posteior parameters --------
- if self.PlotPostDist == True:
-
+ OutputDir = (r'Outputs_Bayes_' + Model.Name + '_' + self.Name)
+ if not os.path.exists(OutputDir): os.makedirs(OutputDir)
+
+ if self.PlotPostDist:
+
import corner
parNames = self.PCEModel.ExpDesign.parNames
figPosterior = corner.corner(self.Posterior_df.to_numpy(), labels=parNames,
@@ -724,11 +724,44 @@ class BayesInference:
figPosterior.set_size_inches((24,16))
- figPosterior.savefig('./'+OutputDir+'/Poisterior_Distribution.svg',
+ if self.emulator:
+ plotname = '/Poisterior_Dist_'+Model.Name+'_emulator'
+ else:
+ plotname = '/Poisterior_Dist_'+Model.Name
+
+ figPosterior.savefig('./'+OutputDir+ plotname + '.svg',
bbox_inches='tight')
+ # -------- Plot logBME dist --------
+ if self.BootstrapItrNr != 1:
+ # Computing the TOM performance
+ x_values = np.linspace(np.min(self.BME), np.max(self.BME), 1000)
+ dof = ObservationData.shape[0]
+ self.logTOMBME = stats.chi2(dof).logpdf(-1*x_values)#pdf(x_values)
+
+ fig, ax = plt.subplots()
+ sns.kdeplot(self.logTOMBME, ax=ax, color="green", shade=True)
+ sns.kdeplot(self.BME, ax=ax, color="blue", shade=True,label='Model BME')
+
+ ax.set_xlabel('log$_{10}$(BME)')
+ ax.set_ylabel('Probability density')
+
+ from matplotlib.patches import Patch
+ legend_elements = [Patch(facecolor='green', edgecolor='green',label='TOM BME'),
+ Patch(facecolor='blue', edgecolor='blue',label='Model BME')]
+ ax.legend(handles=legend_elements,fontsize=14)
+
+ if self.emulator:
+ plotname = '/BME_hist_'+Model.Name+'_emulator'
+ else:
+ plotname = '/BME_hist_'+Model.Name
+
+ plt.savefig('./'+OutputDir+plotname+'.svg', bbox_inches='tight')
+ plt.show()
+ plt.close()
+
# -------- Posteior perdictives --------
- if self.PlotPostPred == True:
+ if self.PlotPostPred:
# sns.set_context("paper")
# Plot the posterior predictive
for Out_idx, OutputName in enumerate(Model.Output.Names):
@@ -737,72 +770,132 @@ class BayesInference:
x_key = list(MeasuredData)[0]
# --- Read prior and posterior predictive ---
- # Prior
- priorPred_df = pd.read_hdf(OutputDir+'/'+"priorPredictive.h5" ,'/data/'+OutputName)
- tags_post = ['prior'] * len(priorPred_df)
- priorPred_df.insert(len(priorPred_df.columns), "Tags", tags_post, True)
+ if self.SamplingMethod == 'rejection':
+ # Prior
+ priorPred_df = pd.read_hdf(OutputDir+'/'+"priorPredictive.h5" ,'/data/'+OutputName)
+ tags_post = ['prior'] * len(priorPred_df)
+ priorPred_df.insert(len(priorPred_df.columns), "Tags", tags_post, True)
- # Posterir
- postPred_df = pd.read_hdf(OutputDir+'/'+"postPredictive.h5" ,'/data/'+OutputName)
- tags_post = ['posterior'] * len(postPred_df)
- postPred_df.insert(len(postPred_df.columns), "Tags", tags_post, True)
+ # Posterior
+ postPred_df = pd.read_hdf(OutputDir+'/'+"postPredictive.h5" ,'/data/'+OutputName)
+ tags_post = ['posterior'] * len(postPred_df)
+ postPred_df.insert(len(postPred_df.columns), "Tags", tags_post, True)
+
+ # Concatenate two dataframes based on x_values
+ frames = [priorPred_df,postPred_df]
+ AllPred_df = pd.concat(frames)
- # Concatenate two dataframes based on x_values
- frames = [priorPred_df,postPred_df]
- AllPred_df = pd.concat(frames)
+ # --- Plot posterior predictive ---
+ sns.violinplot(x_key, y=OutputName, hue="Tags", legend=False,
+ ax=ax, split=True, inner=None, data=AllPred_df, color=".8")
- # --- Plot posterior predictive ---
- sns.violinplot(x_key, y=OutputName, hue="Tags", legend=False,
- ax=ax, split=True, inner=None, data=AllPred_df, color=".8")
+
+ # --- Plot MAP (PCEModel) ---
+ # Find the x,y coordinates for each point
+ x_coords = np.arange(np.unique(AllPred_df[x_key].to_numpy()).shape[0])
+
+ sns.pointplot(x=x_coords, y=MAP_PCEModel[OutputName][0], color='lime',markers='+',
+ linestyles='', capsize=16, ax=ax)
+
+ ax.errorbar(x_coords, MAP_PCEModel[OutputName][0],
+ yerr=1.96*MAP_PCEModelstd[OutputName][0],
+ ecolor='lime', fmt=' ', zorder=-1)
+
+ # --- Plot MAP (origModel) ---
+ sns.pointplot(x=x_coords, y=MAP_origModel[OutputName][0], color='r', markers='+',
+ linestyles='' , capsize=14, ax=ax)
+
+ # --- Plot Data ---
+ ObservationData = MeasuredData
+ first_header = list(ObservationData)[0]
+ ObservationData = ObservationData.round({first_header: 6})
+ sns.pointplot(x=first_header, y=OutputName, color='g', markers='x',
+ linestyles='', capsize=16, data=ObservationData, ax=ax)
+
+ ax.errorbar(x_coords, ObservationData[OutputName].to_numpy(),
+ yerr=1.96*self.MeasurementError[OutputName].to_numpy(),
+ ecolor='g', fmt=' ', zorder=-1)
+
+ # Add labels to the legend
+ handles, labels = ax.get_legend_handles_labels()
+ labels.append('point estimate (PCE Model)')
+ labels.append('point estimate (Orig. Model)')
+ labels.append('Data')
+
+ import matplotlib.lines as mlines
+ #red_patch = mpatches.Patch(color='red')
+ data_marker = mlines.Line2D([], [], color='lime', marker='+', linestyle='None',
+ markersize=10)
+ handles.append(data_marker)
+ data_marker = mlines.Line2D([], [], color='r', marker='+', linestyle='None',
+ markersize=10)
+ handles.append(data_marker)
+ data_marker = mlines.Line2D([], [], color='g', marker='x', linestyle='None',
+ markersize=10)
+ handles.append(data_marker)
+
+ # Add legend
+ ax.legend(handles=handles, labels=labels, loc='best',
+ fontsize='large', frameon=True)
- # --- Plot MAP (PCEModel) ---
- # Find the x,y coordinates for each point
- x_coords = np.arange(np.unique(AllPred_df[x_key].to_numpy()).shape[0])
+ else:
+ # Load posterior predictive
+ AllPred_df = pd.read_hdf(OutputDir+'/'+"postPredictive.h5" ,'/data/'+OutputName)
- ax=sns.pointplot(x=x_coords, y=MAP_PCEModel[OutputName][0], color='lime',markers='+',
- linestyles='', capsize=16, ax=ax) #
-
- ax.errorbar(x_coords, MAP_PCEModel[OutputName][0],
- yerr=1.96*MAP_PCEModelstd[OutputName][0],
- ecolor='lime', fmt=' ', zorder=-1)
-
- # --- Plot MAP (origModel) ---
- sns.pointplot(x=x_coords, y=MAP_origModel[OutputName][0], color='r', markers='+',
- linestyles='' , capsize=14, ax=ax)
-
- # --- Plot Data ---
- ObservationData = MeasuredData
- first_header = list(ObservationData)[0]
- ObservationData = ObservationData.round({first_header: 6})
- sns.pointplot(x=first_header, y=OutputName, color='g', markers='x',
- linestyles='', capsize=16, data=ObservationData, ax=ax)
+ # --- Plot posterior predictive ---
+ x_coords = np.unique(AllPred_df[x_key].to_numpy())
+
+ mu = AllPred_df.groupby(x_key).mean()[key].to_numpy()
+ std = AllPred_df.groupby(x_key).std()[key].to_numpy()
+ plt.plot(x_coords, mu, marker='o', color='b', label='Mean Post. Predictive')
+ plt.fill_between(x_coords, mu-1.96*std, mu+1.96*std, color='b', alpha=0.15)
+
+ # --- Plot Data ---
+ ax.plot(x_coords, MeasuredData[OutputName].to_numpy(),'ko',label ='data',markeredgecolor='w')
+
+ # --- Plot ExpDesign ---
+ index = self.PCEModel.index
+ if index is not None:
+ if self.Name == 'Valid':
+ origOutExpDesign = self.PCEModel.ExpDesign.Y[key][:,index:]
+ else:
+ origOutExpDesign = self.PCEModel.ExpDesign.Y[key][:,:index]
+ else:
+ origOutExpDesign = self.PCEModel.ExpDesign.Y[key]
+
+ for output in origOutExpDesign:
+ plt.plot(x_coords, output, color='grey', alpha=0.15)
+
+ # --- Plot MAP (PCEModel) ---
+ ax.plot(x_coords, MAP_PCEModel[OutputName][0],ls='--', marker='o',c = 'lime' ,label ='point estimate (PCE Model)',markeredgecolor='w')
+
+ # --- Plot MAP (origModel) ---
+ ax.plot(x_coords, MAP_origModel[OutputName][0],ls='-', marker='o',c = 'r',label ='point estimate (Orig. Model)',markeredgecolor='w')
+
+
+ # Add labels for axes
+ plt.xlabel('Time [s]')
+ plt.ylabel(key)
+
+ # Add labels to the legend
+ handles, labels = ax.get_legend_handles_labels()
- ax.errorbar(x_coords, ObservationData[OutputName].to_numpy(),
- yerr=1.96*self.MeasurementError[OutputName].to_numpy(),
- ecolor='g', fmt=' ', zorder=-1)
-
- # Add labels to the legend
- handles, labels = ax.get_legend_handles_labels()
- labels.append('point estimate (PCE Model)')
- labels.append('point estimate (Orig. Model)')
- labels.append('Data')
-
- import matplotlib.lines as mlines
- #red_patch = mpatches.Patch(color='red')
- data_marker = mlines.Line2D([], [], color='lime', marker='+', linestyle='None',
- markersize=10)
- handles.append(data_marker)
- data_marker = mlines.Line2D([], [], color='r', marker='+', linestyle='None',
- markersize=10)
- handles.append(data_marker)
- data_marker = mlines.Line2D([], [], color='g', marker='x', linestyle='None',
- markersize=10)
- handles.append(data_marker)
-
- # Add legend
- ax.legend(handles=handles, labels=labels, loc='best',
- fontsize='large', frameon=True)
-
+ import matplotlib.lines as mlines
+ import matplotlib.patches as mpatches
+ patch = mpatches.Patch(color='b',alpha=0.15)
+ handles.insert(1,patch)
+ labels.insert(1,'95 $\%$ CI')
+
+ data_marker = mlines.Line2D([], [], color='grey',alpha=0.15)
+ handles.append(data_marker)
+ labels.append('Orig. Model Responses')
+
+ # Add legend
+ ax.legend(handles=handles, labels=labels, loc='best',
+ frameon=True)
+
+
+
# Save figure in svg format
if self.emulator:
plotname = '/Post_Prior_Perd_'+Model.Name+'_emulator'
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 5646cbd0e..0a5e8daa4 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -47,7 +47,9 @@ class PostProcessing:
#--------------------------------------------------------------------------------------------------------
def PCEMoments(self, aPCE):
- Model = self.PCEModel.ModelObj
+
+ PCEModel = self.PCEModel
+ Model = PCEModel.ModelObj
for Outkey, ValuesDict in aPCE.CoeffsDict.items():
@@ -62,9 +64,10 @@ class PostProcessing:
PCEMean[idx] = coeffs[0]
# Std = sqrt(sum(coeffs[1:]**2))
PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
-
- self.PCEMeans[Outkey] = PCEMean
- self.PCEStd[Outkey] = PCEStd
+
+ PCA = PCEModel.pca[Outkey]
+ self.PCEMeans[Outkey] = PCA.mean_ + np.dot(PCEMean, PCA.components_) # PCEMean
+ self.PCEStd[Outkey] = np.dot(PCEStd, PCA.components_) #PCEStd
self.Reference = Model.read_MCReference()
@@ -136,7 +139,7 @@ class PostProcessing:
else:
Samples = Samples
self.NrofSamples = len(Samples)
-
+
self.PCEOutputs = {}
self.PCEOutputs_std = {}
univ_p_val = PCEModel.univ_basis_vals(Samples)
@@ -149,33 +152,35 @@ class PostProcessing:
idx = 0
for Inkey, InIdxValues in ValuesDict.items():
- if PCEModel.RegMethod == 'GPE':
- gp = PCEModel.clf_poly[Outkey][Inkey]
- y_mean, y_std = gp.predict(Samples, return_std=True)
+ PolynomialDegrees = PCEModel.BasisDict[Outkey][Inkey]
- else:
- PolynomialDegrees = PCEModel.BasisDict[Outkey][Inkey]
-
- clf_poly = PCEModel.clf_poly[Outkey][Inkey]
-
- PSI_Val = PCEModel.PCE_create_Psi(PolynomialDegrees, univ_p_val)
-
- # Perdiction
- 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)
+ clf_poly = PCEModel.clf_poly[Outkey][Inkey]
+
+ PSI_Val = PCEModel.PCE_create_Psi(PolynomialDegrees, univ_p_val)
+ # Perdiction
+ 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 == '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
-
- self.PCEOutputs[Outkey] = PCEOutputs_mean
- self.PCEOutputs_std[Outkey] = PCEOutputs_std
+
+ PCA = PCEModel.pca[Outkey]
+ self.PCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_) #PCEOutputs_mean
+ self.PCEOutputs_std[Outkey] = np.dot(PCEOutputs_std, PCA.components_) #PCEOutputs_std
return self.PCEOutputs, self.PCEOutputs_std
diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index 71aa79dda..c955d0204 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -218,7 +218,7 @@ class RegressionFastARD():
active = np.zeros(n_features , dtype = np.bool)
- if self.start is not None:
+ if self.start is not None and not hasattr(self, 'active_'):
start = self.start
# start from a given start basis vector
proj = XY**2 / XXd
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 8ff76de96..7270c9885 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -67,6 +67,7 @@ class aPCE:
self.P = None
self.DegreeArray = []
self.RegMethod = None
+ self.metaModel = 'PCE'
self.Discrepancy = None
self.ExpDesign = None
@@ -446,7 +447,7 @@ class aPCE:
return Psi
#--------------------------------------------------------------------------------------------------------
- def RS_Builder(self, PSI, ModelOutput, BasisIndices, RegMethod=None):
+ def RS_Builder(self, PSI, ModelOutput, BasisIndices,startBasisIndices=None, RegMethod=None):
"""
==================================================
@@ -500,7 +501,7 @@ class aPCE:
# n_iter=1000, tol= 0.001)
# alpha_1=1e-6, alpha_2=1e-6,
# lambda_1=Lambda, lambda_2=Lambda)
- clf_poly = RegressionFastARD(fit_intercept=False,compute_score=compute_score,
+ clf_poly = RegressionFastARD(start=startBasisIndices,fit_intercept=False,compute_score=compute_score,
tol= 1e-3)
elif RegMethod == 'LARS':
@@ -539,7 +540,7 @@ class aPCE:
coeffs = clf_poly.coef_
# Raise an error, if all coeffs are zero
- #raise NameError('All the coefficients of the regression model are zero!')
+ #raise ValueError('All the coefficients of the regression model are zero!')
@@ -563,7 +564,7 @@ class aPCE:
return PolynomialDegrees_Sparse,PSI_Sparse,coeffs, clf_poly
#--------------------------------------------------------------------------------------------------------
- def adaptiveSparseaPCE(self, CollocationPoints, ModelOutput, varIdx, verbose=False):
+ def adaptiveSparseaPCE(self, CollocationPoints, ModelOutput, varIdx, prevBasisIndices, verbose=False):
NrSamples,NofPa = CollocationPoints.shape
# Initialization
@@ -615,7 +616,13 @@ class aPCE:
Psi = self.PCE_create_Psi(BasisIndices,univ_p_val)
# ---- Calulate the cofficients of aPCE ----------------------
- sparseBasisIndices, sparsePSI, Coeffs, clf_poly = self.RS_Builder(Psi,ModelOutput,BasisIndices)
+ # Generate the starting indices
+ if prevBasisIndices is not None:
+ t = np.isclose(BasisIndices.toarray()[:, None], prevBasisIndices).all(-1)
+ startIdxes = np.where(t.any(0), t.argmax(0), np.nan).astype(int)
+ else:
+ startIdxes = None
+ sparseBasisIndices, sparsePSI, Coeffs, clf_poly = self.RS_Builder(Psi,ModelOutput,BasisIndices)#,startIdxes)
# Calculate and save the score of LOOCV
score, LCerror = self.CorrectedLeaveoutError(sparsePSI, Coeffs, ModelOutput, clf_poly)
@@ -824,7 +831,7 @@ class aPCE:
Q_2 = 1 - CorrectedErrLoo
- return Q_2, LCerror
+ return Q_2, residual #LCerror
#--------------------------------------------------------------------------------------------------------
# Initialization
@@ -856,7 +863,7 @@ class aPCE:
gp.fit(X, y)
# Compute score
- Score = gp.score(X, y)
+ # Score = gp.score(X, y)
# print('\n'+'-'*50)
# print('The estimation of GPE coefficients converged at for output variable %s.'%(varIdx+1))
@@ -885,39 +892,41 @@ class aPCE:
self.BasisDict = self.AutoVivification()
self.ScoresDict = self.AutoVivification()
self.clf_poly = self.AutoVivification()
- # self.gp_poly = self.AutoVivification()
+ self.gp_poly = self.AutoVivification()
+ self.pca = self.AutoVivification()
self.LCerror = self.AutoVivification()
self.allBasisIndices = self.AutoVivification()
# Read observations or MCReference
if Model.MeasurementFile is not None or len(Model.Observations) != 0:
self.Observations = Model.read_Observation()
-
- # self.DegreeArray = np.array(range(self.MinPceDegree,self.MaxPceDegree+1))
- # self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
- # for degIdx, deg in enumerate(self.DegreeArray):
- # for qidx, q in enumerate(self.q):
- # # Generate the polynomial basis indices
- # self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree = deg, q=q)
-
-
- # Generate Basis indices matrix
- maxDeg = self.MaxPceDegree
- nSamples, ndim = self.OptimalCollocationPointsBase.shape
- nitr = nSamples - self.ExpDesign.initNrSamples
- d = nitr if nitr != 0 and self.NofPa > 5 else 1
- M_uptoMax = lambda maxDeg: np.array([math.factorial(ndim+d)//(math.factorial(ndim)*math.factorial(d)) for d in range(1,maxDeg+1)])
- deg = range(1,maxDeg+1)[np.argmin(abs(M_uptoMax(maxDeg)-nSamples*ndim*d))]
- self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
- if deg not in self.DegreeArray:
- self.allBasisIndices = self.AutoVivification()
- self.DegreeArray = np.arange(self.MinPceDegree,deg+1)
+ if self.ExpDesign.Method == 'normal':
+ self.DegreeArray = np.array(range(self.MinPceDegree,self.MaxPceDegree+1))
+ self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
for degIdx, deg in enumerate(self.DegreeArray):
for qidx, q in enumerate(self.q):
# Generate the polynomial basis indices
self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree=deg, q=q)
+ else:
+ # Generate Basis indices matrix
+ maxDeg = self.MaxPceDegree
+ nSamples, ndim = self.OptimalCollocationPointsBase.shape
+ nitr = nSamples - self.ExpDesign.initNrSamples
+ d = nitr if nitr != 0 and self.NofPa > 5 else 1
+ M_uptoMax = lambda maxDeg: np.array([math.factorial(ndim+d)//(math.factorial(ndim)*math.factorial(d)) for d in range(1,maxDeg+1)])
+ deg = range(1,maxDeg+1)[np.argmin(abs(M_uptoMax(maxDeg)-nSamples*ndim*d))]
+ self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
+
+ if deg not in self.DegreeArray:
+ self.allBasisIndices = self.AutoVivification()
+ self.DegreeArray = np.arange(self.MinPceDegree,deg+1)
+ for degIdx, deg in enumerate(self.DegreeArray):
+ for qidx, q in enumerate(self.q):
+ # Generate the polynomial basis indices
+ self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree=deg, q=q)
+
# Evaluate the univariate polynomials on ExpDesign
if self.RegMethod != 'GPE':
@@ -927,48 +936,71 @@ class aPCE:
if 'x_values' in OutputDict: del OutputDict['x_values']
# Loop through data points and fit the surrogate
- for key , OutputMatrix in OutputDict.items():
-
- if not OptDesignFlag: print(">>>> Calculating the PCE coefficients of %s with %s<<<<<<"%(key, self.RegMethod))
- for idx in range(OutputMatrix.shape[1]):
+ for key , Output in OutputDict.items():
+
+ # TODO: Transform via Principal Component Analysis
+ from sklearn.decomposition import PCA as sklearnPCA
+ n_samples, n_features = Output.shape
+ covar_matrix = sklearnPCA(n_components=None)
+ covar_matrix.fit(Output)
+ var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_, decimals=5)*100)
+ try:
+ selected_n_components = np.where(var>=99.999)[0][0] + 1
+ except:
+ selected_n_components = min(n_samples, n_features)
- # Only update the coeffs in the searching alg. of SeqExpDesign
- if OptDesignFlag:
- BasisIndices = self.BasisDict[key]["y_"+str(idx+1)]
- self.CoeffsDict[key]["y_"+str(idx+1)] , self.BasisDict[key]["y_"+str(idx+1)], self.ScoresDict[key]["y_"+str(idx+1)], self.clf_poly[key]["y_"+str(idx+1)], self.LCerror[key]["y_"+str(idx+1)] = self.CoeffsUpdateaPCE(CollocationPoints, OutputMatrix[:,idx], BasisIndices)
-
- else:
- # TODO: Basis expansion
- # if self.ExpDesignFlag == 'sequential':
- # self.allBasisIndices = self.AutoVivification()
- # BasisIndices = self.BasisDict[key]["y_"+str(idx+1)]
- # lastmaxdeg = np.max(BasisIndices)#np.max(np.sum(BasisIndices,axis=1))
- # T = 1
- # for t in np.arange(1,T+1):
- # q = 1.0
- # for indices in BasisIndices:
- # for j, index in enumerate(indices):
- # lambda_n = indices.copy()
- # lambda_n[j] += 1 #
- # if sum(lambda_n)>lastmaxdeg and sum(lambda_n)<=lastmaxdeg+t \
- # and not lambda_n.tolist() in BasisIndices.tolist()\
- # and all(np.array(lambda_n)<=self.MaxPceDegree):
- # # if not lambda_n.tolist() in BasisIndices.tolist()\
- # # and all(np.array(lambda_n)<=self.MaxPceDegree):
- # BasisIndices = np.vstack((BasisIndices,lambda_n))
-
- # self.allBasisIndices[str(t)][str(q)] = BasisIndices
+ nComponents = min(n_samples, n_features, selected_n_components)
+
+ self.pca[key] = sklearnPCA(n_components=nComponents)
+ OutputMatrix = self.pca[key].fit_transform(Output)
+
+ if self.RegMethod == 'GPE': # TODO: Use multi response GPE
+ self.gp_poly[key] = self.GaussianProcessEmulator(CollocationPoints, OutputMatrix)
+
+ else:
+ if not OptDesignFlag: print(">>>> Calculating the PCE coefficients of %s with %s<<<<<<"%(key, self.RegMethod))
+ for idx in range(OutputMatrix.shape[1]):
- self.CoeffsDict[key]["y_"+str(idx+1)] , self.BasisDict[key]["y_"+str(idx+1)], self.ScoresDict[key]["y_"+str(idx+1)], self.clf_poly[key]["y_"+str(idx+1)], self.LCerror[key]["y_"+str(idx+1)] = self.adaptiveSparseaPCE(CollocationPoints,
- OutputMatrix[:,idx],idx,verbose=True)
+ # Only update the coeffs in the searching alg. of SeqExpDesign
+ if OptDesignFlag:
+ BasisIndices = self.BasisDict[key]["y_"+str(idx+1)]
+ self.CoeffsDict[key]["y_"+str(idx+1)] , self.BasisDict[key]["y_"+str(idx+1)], self.ScoresDict[key]["y_"+str(idx+1)], self.clf_poly[key]["y_"+str(idx+1)], self.LCerror[key]["y_"+str(idx+1)] = self.CoeffsUpdateaPCE(CollocationPoints, OutputMatrix[:,idx], BasisIndices)
+
+ else:
+ # TODO: Basis expansion
+ # if self.ExpDesignFlag == 'sequential':
+ # self.allBasisIndices = self.AutoVivification()
+ # BasisIndices = self.BasisDict[key]["y_"+str(idx+1)]
+ # lastmaxdeg = np.max(BasisIndices)#np.max(np.sum(BasisIndices,axis=1))
+ # T = 1
+ # for t in np.arange(1,T+1):
+ # q = 1.0
+ # for indices in BasisIndices:
+ # for j, index in enumerate(indices):
+ # lambda_n = indices.copy()
+ # lambda_n[j] += 1 #
+ # if sum(lambda_n)>lastmaxdeg and sum(lambda_n)<=lastmaxdeg+t \
+ # and not lambda_n.tolist() in BasisIndices.tolist()\
+ # and all(np.array(lambda_n)<=self.MaxPceDegree):
+ # # if not lambda_n.tolist() in BasisIndices.tolist()\
+ # # and all(np.array(lambda_n)<=self.MaxPceDegree):
+ # BasisIndices = np.vstack((BasisIndices,lambda_n))
+
+ # self.allBasisIndices[str(t)][str(q)] = BasisIndices
+
+ prevBasisIndices = self.BasisDict[key]["y_"+str(idx)] if idx != 0 else None
+
+
+ self.CoeffsDict[key]["y_"+str(idx+1)] , self.BasisDict[key]["y_"+str(idx+1)], self.ScoresDict[key]["y_"+str(idx+1)], self.clf_poly[key]["y_"+str(idx+1)], self.LCerror[key]["y_"+str(idx+1)] = self.adaptiveSparseaPCE(CollocationPoints,OutputMatrix[:,idx],
+ idx,prevBasisIndices,verbose=True)
- # TODO: PCEKriging (Gaussian Process Emulator)
- # if self.RegMethod == 'PCEKriging':
- # self.gp_poly[key]["y_"+str(idx+1)] = self.GaussianProcessEmulator(CollocationPoints, self.LCerror[key]["y_"+str(idx+1)])
-
+ # TODO: PCEKriging (PCE with Gaussian Process Emulator)
+ if self.metaModel == 'PCEKriging':
+ self.gp_poly[key]["y_"+str(idx+1)] = self.GaussianProcessEmulator(CollocationPoints, self.LCerror[key]["y_"+str(idx+1)])
+
return self
#--------------------------------------------------------------------------------------------------------
@@ -1161,7 +1193,8 @@ class aPCE:
# Enlarge sample size if it doesn't fulfill the criteria
if ((ESS > MCsize) or (ESS < 1)):
- MCsize = MCsize * 10
+ MCsize *= 10
+ ESS = 0
# Rejection Step
# Random numbers between 0 and 1
@@ -1170,7 +1203,6 @@ class aPCE:
# Reject the poorly performed prior
accepted = (Likelihoods/np.max(Likelihoods)) >= unif
-
# Prior-based estimation of BME
logBME = np.log(np.nanmean(Likelihoods))
@@ -1289,7 +1321,8 @@ class aPCE:
# Enlarge sample size if it doesn't fulfill the criteria
if ((ESS > MCsize) or (ESS < 1)):
- MCsize = MCsize * 10
+ MCsize *= 10
+ ESS = 0
# Rejection Step
# Random numbers between 0 and 1
@@ -1942,11 +1975,10 @@ class aPCE:
#--------------------------------------------------------------------------------------------------------
def eval_PCEmodel(self,Samples):
- meanPCEOutputs = {}
- stdPCEOutputs = {}
univ_p_val = self.univ_basis_vals(Samples)
-
+ meanPCEOutputs = {}
+ stdPCEOutputs = {}
for Outkey, ValuesDict in self.CoeffsDict.items():
PCEOutputs_mean = np.zeros((len(Samples), len(ValuesDict)))
@@ -1955,32 +1987,31 @@ class aPCE:
for Inkey, InIdxValues in ValuesDict.items():
idx = int(Inkey.split('_')[1]) - 1
- if self.RegMethod == 'GPE':
- gp = self.clf_poly[Outkey][Inkey]
- PCEOutputs_mean[:, idx], PCEOutputs_std[:, idx] = gp.predict(Samples, return_std=True)
+ PolynomialDegrees = self.BasisDict[Outkey][Inkey]
+ PSI_Val = self.PCE_create_Psi(PolynomialDegrees, univ_p_val)
- else:
-
- 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)
- # Perdiction
- try: # with error bar
- # gp_poly = self.gp_poly[Outkey][Inkey]
- # gp_mean,PCEOutputs_std[:, idx] = gp_poly.predict(Samples, return_std=True)
- clf_poly = self.clf_poly[Outkey][Inkey]
- PCEOutputs_mean[:, idx], PCEOutputs_std[:, idx] = clf_poly.predict(PSI_Val, return_std=True)
- # PCEOutputs_mean[:, idx] = clf_poly.predict(PSI_Val)
- # PCEOutputs_mean[:, idx] += gp_mean
- # PCEOutputs_mean[:, idx] = y_mean
- # PCEOutputs_std[:, idx] = y_std
-
- except: # without error bar
- coeffs = self.CoeffsDict[Outkey][Inkey]
- PCEOutputs_mean[:, idx] = np.dot(PSI_Val, coeffs)
+ 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 == 'PCEKriging':
+ 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
+
+ PCEOutputs_mean[:, idx] = y_mean
+ PCEOutputs_std[:, idx] = y_std
- meanPCEOutputs[Outkey] = PCEOutputs_mean
- stdPCEOutputs[Outkey] = PCEOutputs_std
+ PCA = self.pca[Outkey]
+ meanPCEOutputs[Outkey] = PCA.mean_ + np.dot(PCEOutputs_mean,PCA.components_) #PCEOutputs_mean
+ stdPCEOutputs[Outkey] = np.dot(PCEOutputs_std, PCA.components_) #PCEOutputs_std
return meanPCEOutputs, stdPCEOutputs
@@ -2004,7 +2035,7 @@ class aPCE:
covMatrix = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
np.fill_diagonal(covMatrix, Sigma2s)
- # Add the std of the PCE is emulator is chosen.
+ # Add the std of the PCE.
covMatrix_PCE = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
stdPCE = np.empty((SampleSize,0))
for key in ModelOutputNames:
@@ -2012,8 +2043,10 @@ class aPCE:
# Expected value of variance (Assump: i.i.d stds)
varPCE = np.mean(stdPCE**2, axis=0)
-
+
np.fill_diagonal(covMatrix_PCE, varPCE)
+
+ # Aggregate the cov matrices
covMatrix = covMatrix + covMatrix_PCE
# Compute likelihood
@@ -2377,6 +2410,10 @@ class aPCE:
# Std = sqrt(sum(coeffs[1:]**2))
PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
+ PCA = self.pca[Outkey]
+ PCEMean = PCA.mean_ + np.dot(PCEMean, PCA.components_) # PCEMean
+ PCEStd = np.dot(PCEStd, PCA.components_) #PCEStd
+
# Compute the error between mean and std of PCEModel and OrigModel
RMSE_Mean = mean_absolute_error(df_MCReference['mean'], PCEMean)
#np.sqrt(mean_squared_error(df_MCReference['mean'], PCEMean))
@@ -2452,7 +2489,9 @@ class aPCE:
Scores_all, varExpDesignY =[], []
for OutName in OutputName:
Scores_all.append(list(initPCEModel.ScoresDict[OutName].values()))
- varExpDesignY.append(np.var(initPCEModel.ExpDesign.Y[OutName],axis=0))
+ components = self.pca[OutName].transform(initPCEModel.ExpDesign.Y[OutName])
+ varExpDesignY.append(np.var(components,axis=0))
+ # varExpDesignY.append(np.var(initPCEModel.ExpDesign.Y[OutName],axis=0))
Scores = [item for sublist in Scores_all for item in sublist]
weights = [item for sublist in varExpDesignY for item in sublist]
initModifiedLOO = [np.average([1-score for score in Scores], weights=weights)]
@@ -2565,7 +2604,9 @@ class aPCE:
Scores_all, varExpDesignY =[], []
for OutName in OutputName:
Scores_all.append(list(PCEModel.ScoresDict[OutName].values()))
- varExpDesignY.append(np.var(PCEModel.ExpDesign.Y[OutName],axis=0))
+ components = self.pca[OutName].transform(initPCEModel.ExpDesign.Y[OutName])
+ varExpDesignY.append(np.var(components,axis=0))
+ # varExpDesignY.append(np.var(PCEModel.ExpDesign.Y[OutName],axis=0))
Scores = [item for sublist in Scores_all for item in sublist]
weights = [item for sublist in varExpDesignY for item in sublist]
ModifiedLOO = [np.average([1-score for score in Scores],
@@ -2729,7 +2770,6 @@ class aPCE:
for key in list(OutputDict.keys())[1:]:
OutputMatrix = OutputDict[key]
- print("--- Slicing of %s ---"%key)
for idx in range(OutputMatrix.shape[1]):
newPCEModel.CoeffsDict[key]["y_"+str(index+idx+1)] = PCEModel.CoeffsDict[key]["y_"+str(index+idx+1)]
newPCEModel.BasisDict[key]["y_"+str(index+idx+1)] = PCEModel.BasisDict[key]["y_"+str(index+idx+1)]
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
index 29a27c08d..c47ea65c3 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalytFuncValid_Test.py
@@ -48,7 +48,7 @@ if __name__ == "__main__":
# For Bayesian inversion synthetic data with X=[0,0]
Model.Observations['Time [s]'] = np.arange(0, 10, 1.) / 5
- Model.Observations['Z'] = np.array([[2]*5]).reshape(-1,1)
+ Model.Observations['Z'] = np.array([[2]*5])
# For Checking with the MonteCarlo refrence
# Model.MCReference['Time [s]'] = np.arange(0, 10, 1.) / 9
@@ -115,18 +115,17 @@ if __name__ == "__main__":
NrofInitSamples = 25
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
+ MetaModelOpts.ExpDesign.Method = 'normal' # 1) normal 2) sequential
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.MaxNSamples = 50
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-15
# Defining the measurement error, if it's known a priori
- obsData = pd.DataFrame(Model.Observations, columns=Model.Output.Names)**2
-
+ obsData = pd.DataFrame(Model.Observations, columns=Model.Output.Names)
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = obsData
+ DiscrepancyOpts.Parameters = obsData ** 2
MetaModelOpts.Discrepancy = DiscrepancyOpts
# Plot the posterior snapshots for SeqDesign
@@ -143,7 +142,7 @@ if __name__ == "__main__":
#MetaModelOpts.ExpDesign.nReprications = 2
# -------- Exploration ------
#1)'Voronoi' 2)'MC' 3)'LHS' 4)'dual annealing'
- MetaModelOpts.ExpDesign.ExploreMethod = 'NewMethod'
+ MetaModelOpts.ExpDesign.ExploreMethod = 'Voronoi'
# Use when 'dual annealing' chosen
MetaModelOpts.ExpDesign.MaxFunItr = 200
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index b6d52faa7..9a82e0aa5 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -92,8 +92,8 @@ if __name__ == "__main__":
# The degree with the lowest Leave-One-Out cross-validation (LOO)
# error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
- # MetaModelOpts.MinPceDegree = 1 #12
- MetaModelOpts.MaxPceDegree = 15 #12
+ MetaModelOpts.MinPceDegree = 1 #12
+ MetaModelOpts.MaxPceDegree = 12 #12
# q-quasi-norm 0<q<1 (default=1)
MetaModelOpts.q = 1.0 if ndim<5 else 0.5
@@ -103,7 +103,7 @@ if __name__ == "__main__":
# 1)OLS: Ordinary Least Square 2)BRR: Bayesian Ridge Regression
# 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression
# 5)SGDR: Stochastic gradient descent learning
-
+ # MetaModelOpts.metaModel = 'PCEKriging'
MetaModelOpts.RegMethod = 'ARD'
# Print summary of the regression results
@@ -116,7 +116,7 @@ if __name__ == "__main__":
# One-shot (normal) or Sequential Adaptive (sequential) Design
MetaModelOpts.ExpDesign.Method = 'sequential'
- NrofInitSamples = 100#75
+ NrofInitSamples = 100 #75
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
# Sampling methods
@@ -146,7 +146,7 @@ if __name__ == "__main__":
# ------- Sequential Design configuration --------
# ------------------------------------------------
# 1) 'None' 2) 'equal' 3)'epsilon-decreasing' 4) 'adaptive'
- MetaModelOpts.ExpDesign.TradeOffScheme = 'None'
+ MetaModelOpts.ExpDesign.TradeOffScheme = 'epsilon-decreasing'
# MetaModelOpts.ExpDesign.nReprications = 2
# -------- Exploration ------
#1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'LOOCV' 5)'dual annealing'
@@ -167,7 +167,7 @@ if __name__ == "__main__":
# 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 = 'AIC'
+ MetaModelOpts.ExpDesign.UtilityFunction = 'DKL'
# VarBasedOptDesign -> when data is not available
# Only with Vornoi >>> 1)Entropy 2)EIGF, 3)LOOCV
@@ -223,6 +223,7 @@ if __name__ == "__main__":
#======== Bayesian inference with Emulator ==========
#=====================================================
BayesOpts = BayesInference(PCEModel)
+ # BayesOpts.SamplingMethod = 'MCMC'
BayesOpts.emulator = True
# Bootstrap
@@ -242,7 +243,7 @@ if __name__ == "__main__":
# (Option B)
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = obsData **2
+ DiscrepancyOpts.Parameters = obsData**2
BayesOpts.Discrepancy = DiscrepancyOpts
@@ -259,7 +260,8 @@ if __name__ == "__main__":
BayesOptsModel.NrofSamples = 100000
# Evaluation of forward model predictions for the samples used for PCEModel
- BayesOptsModel.Samples = Bayes_PCE.Samples
+ # BayesOptsModel.Samples = Bayes_PCE.Samples
+ BayesOptsModel.SamplingMethod = 'MCMC'
# Bootstrap for BME calulations
BayesOptsModel.Bootstrap = True
--
GitLab