diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 9061ee1beb5b12b1509bb614ced7713a228b9766..9b8de34b44874345d5acb9e5542f936051263aad 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -259,7 +259,6 @@ class BayesInference:
# Add the std of the PCE is chosen as emulator.
if self.emulator:
stdPCE = std[out] if std is not None else np.mean(self._stdPCEPriorPred[out], axis=0)
- # stdPCE = np.mean(self._stdPCEPriorPred[out], axis=0)
# Expected value of variance (Assump: i.i.d stds)
covMatrix += np.diag(stdPCE**3)
@@ -648,6 +647,7 @@ class BayesInference:
if self.emulator:
stdPCE = PCEModel.RMSE[var] if PCEModel.RMSE is not None else PosteriorPred_std[var][i]
+ stdPCE = PosteriorPred_std[var][i]
# Expected value of variance (Assump: i.i.d stds)
cov += np.diag(stdPCE**3)
@@ -851,7 +851,7 @@ class BayesInference:
# Train a GPR meta-model using MAP
self.errorMetaModel = PCEModel.create_ModelError(self.BiasInputs,y_MAP,
- Name=self.Name)
+ Name=self.Name)
# -----------------------------------------------------
# ----- Loop over the perturbed observation data ------
# -----------------------------------------------------
@@ -869,7 +869,7 @@ class BayesInference:
name=self.Name)
# TODO: Correct the predictions with Model discrepancy
if hasattr(self, 'errorModel') and self.errorModel:
- self.__meanPCEPriorPred, PosteriorPred_std = self.errorMetaModel.eval_errormodel(self.BiasInputs,
+ self.__meanPCEPriorPred, self._stdPCEPriorPred = self.errorMetaModel.eval_errormodel(self.BiasInputs,
self.__meanPCEPriorPred)
# Surrogate model's error using RMSE of test data
diff --git a/BayesValidRox/BayesInference/MCMC.py b/BayesValidRox/BayesInference/MCMC.py
index c04871b3e1aa214c1199f43d415a15195fc6c158..77efc6954181046d4ffe2b0d7c4c28207141079a 100755
--- a/BayesValidRox/BayesInference/MCMC.py
+++ b/BayesValidRox/BayesInference/MCMC.py
@@ -451,7 +451,7 @@ class MCMC():
self.counter += 1
if hasattr(self, 'errorMetaModel') and BayesObj.errorModel:
- meanPred, _ = self.errorMetaModel.eval_errormodel(BayesObj.BiasInputs,meanPred)
+ meanPred, stdPred = self.errorMetaModel.eval_errormodel(BayesObj.BiasInputs,meanPred)
return meanPred, stdPred
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 337440c7147e9246f1a920e28537bbf2cb5c1d95..a014e709a548bee3c8323eac8779e70a7d4eff02 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -16,6 +16,7 @@ import h5py
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnchoredText
+from sklearn.preprocessing import MinMaxScaler
SIZE = 30
plt.rc('figure', figsize = (24, 16))
@@ -41,6 +42,11 @@ from scipy import stats
from sklearn.linear_model import LinearRegression, ARDRegression
from sklearn import linear_model
from sklearn.metrics import r2_score, mean_squared_error
+from sklearn.gaussian_process import GaussianProcessRegressor
+import sklearn.gaussian_process.kernels as kernels
+# import (RBF, Matern, RationalQuadratic,
+# ExpSineSquared, DotProduct,
+# ConstantKernel)
import os
import sys
import multiprocessing
@@ -51,6 +57,7 @@ from scipy.spatial import distance
import scipy.linalg as spla
from numpy.polynomial.polynomial import polyval
import chaospy
+from joblib import Parallel,delayed
from BayesInference.Discrepancy import Discrepancy
from .ExpDesigns import ExpDesigns
@@ -718,23 +725,18 @@ class Metamodel:
#--------------------------------------------------------------------------------------------------------
def GaussianProcessEmulator(self, X, y, nugTerm=None,autoSelect=False, varIdx=None):
- from sklearn.gaussian_process import GaussianProcessRegressor
- from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic,
- ExpSineSquared, DotProduct,
- ConstantKernel)
-
nugTerm = nugTerm if nugTerm else np.var(y)
- kernels = [nugTerm * RBF(length_scale=1.0, length_scale_bounds=(1e-15, 1e5)),
- nugTerm * RationalQuadratic(length_scale=0.2, alpha=1.0),
- # 1.0 * ExpSineSquared(length_scale=1.0, length_scale_bounds=(1e-05, 1e05)),
- nugTerm * Matern(length_scale=1.0, length_scale_bounds=(1e-1, 1e1),
- nu=1.5)]
+ Kernels = [ nugTerm * kernels.RBF(length_scale=1.0, length_scale_bounds=(1e-25, 1e15)),
+ nugTerm * kernels.RationalQuadratic(length_scale=0.2, alpha=1.0),
+ # 1.0 * kernels.ExpSineSquared(length_scale=1.0, length_scale_bounds=(1e-05, 1e05)),
+ nugTerm * kernels.Matern(length_scale=1.0, length_scale_bounds=(1e-1, 1e1),
+ nu=2.5)]
if autoSelect:# Automatic selection of the kernel
gp = {}
BME = []
- for i, kernel in enumerate(kernels):
+ for i, kernel in enumerate(Kernels):
gp[i] = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=2,
normalize_y=True)
@@ -746,7 +748,7 @@ class Metamodel:
gp = gp[np.argmax(BME)]
else:
- gp = GaussianProcessRegressor(kernel=kernels[0], n_restarts_optimizer=5,
+ gp = GaussianProcessRegressor(kernel=Kernels[0], n_restarts_optimizer=3,
normalize_y=False)
gp.fit(X, y)
@@ -839,8 +841,6 @@ class Metamodel:
# Parallel fit regression
- from joblib import Parallel,delayed
-
if self.metaModel.lower() == 'gpe':
outDict = Parallel(n_jobs=-1, prefer='threads')(
delayed(self.GaussianProcessEmulator)(CollocationPoints,OutputMatrix[:,idx])
@@ -894,18 +894,8 @@ class Metamodel:
Model = self.ModelObj
outputNames = Model.Output.Names
self.errorRegMethod = 'GPE'
- self.errorCoeffsDict= self.AutoVivification()
- self.errorBasisDict = self.AutoVivification()
- self.errorScoresDict = self.AutoVivification()
self.errorclf_poly = self.AutoVivification()
- self.errorLCerror = self.AutoVivification()
self.errorScale = self.AutoVivification()
- # Original EDY
- EDY = self.ExpDesign.Y
-
- # Prepare the input matrix
- from sklearn.preprocessing import MinMaxScaler
- # self.errorScale = {}
# Read data
if Name.lower() == 'calib':
@@ -936,7 +926,7 @@ class Metamodel:
def eval_errormodel(self,X,y_pred):
meanPCEOutputs = {}
stdPCEOutputs = {}
-
+
for Outkey, ValuesDict in self.errorclf_poly.items():
PCEOutputs_mean = np.zeros_like(y_pred[Outkey])
@@ -951,7 +941,7 @@ class Metamodel:
for j,pred in enumerate(y_pred[Outkey]):
BiasInputs = X[Outkey]#np.hstack((X[Outkey],pred.reshape(-1,1)))
Samples_S = scaler.transform(BiasInputs)
-
+
PCEOutputs_mean[j], PCEOutputs_std[j] = gp.predict(Samples_S, return_std=True)
PCEOutputs_mean[j] += pred