diff --git a/src/bayesvalidrox/surrogate_models/meta_model_engine.py b/src/bayesvalidrox/surrogate_models/meta_model_engine.py
index 08d91ecf77ff20aaaf88384604d75e22b31f8cc7..3985da8a3d4ff124428c9f8710b4d029ca03716c 100644
--- a/src/bayesvalidrox/surrogate_models/meta_model_engine.py
+++ b/src/bayesvalidrox/surrogate_models/meta_model_engine.py
@@ -486,7 +486,7 @@ class MetaModelEngine():
return -1 * score # -1 is for minimization instead of maximization
# -------------------------------------------------------------------------
- def util_BayesianActiveDesign(self, X_can, sigma2Dict, var='DKL'):
+ def util_BayesianActiveDesign(self, y_hat, std, sigma2Dict, var='DKL'):
"""
Computes scores based on Bayesian active design criterion (var).
@@ -511,98 +511,35 @@ class MetaModelEngine():
"""
- # Evaluate the PCE metamodels at that location ???
- Y_mean_can, Y_std_can = self.MetaModel.eval_metamodel(
- samples=np.array([X_can])
- )
-
# Get the data
obs_data = self.observations
n_obs = self.Model.n_obs
- # TODO: Analytical DKL
+ mc_size = 10000
+
# Sample a distribution for a normal dist
# with Y_mean_can as the mean and Y_std_can as std.
-
- # priorMean, priorSigma2, Obs = np.empty((0)),np.empty((0)),np.empty((0))
-
- # for key in list(Y_mean_can):
- # # concatenate the measurement error
- # Obs = np.hstack((Obs,ObservationData[key]))
-
- # # concatenate the mean and variance of prior predictive
- # means, stds = Y_mean_can[key][0], Y_std_can[key][0]
- # priorMean = np.hstack((priorSigma2,means))
- # priorSigma2 = np.hstack((priorSigma2,stds**2))
-
- # # Covariance Matrix of prior
- # covPrior = np.zeros((priorSigma2.shape[0], priorSigma2.shape[0]), float)
- # np.fill_diagonal(covPrior, priorSigma2)
-
- # # Covariance Matrix of Likelihood
- # covLikelihood = np.zeros((sigma2Dict.shape[0], sigma2Dict.shape[0]), float)
- # np.fill_diagonal(covLikelihood, sigma2Dict)
-
- # # Calculate moments of the posterior (Analytical derivation)
- # n = priorSigma2.shape[0]
- # covPost = np.dot(np.dot(covPrior,np.linalg.inv(covPrior+(covLikelihood/n))),covLikelihood/n)
-
- # meanPost = np.dot(np.dot(covPrior,np.linalg.inv(covPrior+(covLikelihood/n))) , Obs) + \
- # np.dot(np.dot(covPrior,np.linalg.inv(covPrior+(covLikelihood/n))),
- # priorMean/n)
- # # Compute DKL from prior to posterior
- # term1 = np.trace(np.dot(np.linalg.inv(covPrior),covPost))
- # deltaMean = priorMean-meanPost
- # term2 = np.dot(np.dot(deltaMean,np.linalg.inv(covPrior)),deltaMean[:,None])
- # term3 = np.log(np.linalg.det(covPrior)/np.linalg.det(covPost))
- # DKL = 0.5 * (term1 + term2 - n + term3)[0]
-
- # ---------- Inner MC simulation for computing Utility Value ----------
- # Estimation of the integral via Monte Varlo integration
- MCsize = 20000
- ESS = 0
-
- while ((ESS > MCsize) or (ESS < 1)):
-
- # Sample a distribution for a normal dist
- # with Y_mean_can as the mean and Y_std_can as std.
- Y_MC, std_MC = {}, {}
- logPriorLikelihoods = np.zeros((MCsize))
- for key in list(Y_mean_can):
- means, stds = Y_mean_can[key][0], Y_std_can[key][0]
- # cov = np.zeros((means.shape[0], means.shape[0]), float)
- # np.fill_diagonal(cov, stds**2)
-
- Y_MC[key] = np.zeros((MCsize, n_obs))
- logsamples = np.zeros((MCsize, n_obs))
- for i in range(n_obs):
- NormalDensity = stats.norm(means[i], stds[i])
- Y_MC[key][:, i] = NormalDensity.rvs(MCsize)
- logsamples[:, i] = NormalDensity.logpdf(Y_MC[key][:, i])
-
- logPriorLikelihoods = np.sum(logsamples, axis=1)
- std_MC[key] = np.zeros((MCsize, means.shape[0]))
-
- # Likelihood computation (Comparison of data and simulation
- # results via PCE with candidate design)
- likelihoods = self.__normpdf(Y_MC, std_MC, obs_data, sigma2Dict)
-
- # Check the Effective Sample Size (1<ESS<MCsize)
- ESS = 1 / np.sum(np.square(likelihoods/np.nansum(likelihoods)))
-
- # Enlarge sample size if it doesn't fulfill the criteria
- if ((ESS > MCsize) or (ESS < 1)):
- MCsize *= 10
- ESS = 0
+ Y_MC, std_MC = {}, {}
+ logPriorLikelihoods = np.zeros((mc_size))
+ for key in list(y_hat):
+ cov = np.diag(std[key]**2)
+ rv = stats.multivariate_normal(mean=y_hat[key], cov=cov)
+ Y_MC[key] = rv.rvs(size=mc_size)
+ logPriorLikelihoods += rv.logpdf(Y_MC[key])
+ std_MC[key] = np.zeros((mc_size, y_hat[key].shape[0]))
+
+ # Likelihood computation (Comparison of data and simulation
+ # results via PCE with candidate design)
+ likelihoods = self.__normpdf(Y_MC, std_MC, obs_data, sigma2Dict)
# Rejection Step
# Random numbers between 0 and 1
- unif = np.random.rand(1, MCsize)[0]
+ unif = np.random.rand(1, mc_size)[0]
# Reject the poorly performed prior
accepted = (likelihoods/np.max(likelihoods)) >= unif
# Prior-based estimation of BME
- logBME = np.log(np.nanmean(likelihoods))
+ logBME = np.log(np.nanmean(likelihoods), dtype=np.float128)
# Posterior-based expectation of likelihoods
postLikelihoods = likelihoods[accepted]
@@ -626,7 +563,7 @@ class MetaModelEngine():
# Marginal log likelihood
elif var == 'BME':
- U_J_d = logBME
+ U_J_d = np.nanmean(likelihoods)
# Entropy-based information gain
elif var == 'infEntropy':
@@ -656,7 +593,7 @@ class MetaModelEngine():
# D_theta_bar = np.mean(-2 * Likelihoods)
N_star_p = 0.5 * np.var(np.log(likelihoods[likelihoods != 0]))
Likelihoods_theta_mean = self.__normpdf(
- Y_mean_can, Y_std_can, obs_data, sigma2Dict
+ y_hat, std, obs_data, sigma2Dict
)
DIC = -2 * np.log(Likelihoods_theta_mean) + 2 * N_star_p
@@ -821,8 +758,7 @@ class MetaModelEngine():
"""
# To avoid changes ub original aPCE object
- Model = self.Model
- MetaModel = deepcopy(self.MetaModel)
+ MetaModel = self.MetaModel
out_names = MetaModel.ModelObj.Output.names
if X_can.ndim == 1:
X_can = X_can.reshape(1, -1)
@@ -838,25 +774,29 @@ class MetaModelEngine():
# Evaluate the PCE metamodels at that location ???
Y_PC_can, Y_std_can = MetaModel.eval_metamodel(samples=X_can)
-
- # Add all suggestion as new ExpDesign
+ PCE_Model_can = deepcopy(MetaModel)
+ # Add the candidate to the ExpDesign
NewExpDesignX = np.vstack((oldExpDesignX, X_can))
NewExpDesignY = {}
for key in oldExpDesignY.keys():
- try:
- NewExpDesignY[key] = np.vstack((oldExpDesignY[key],
- Y_PC_can[key]))
- except:
- NewExpDesignY[key] = oldExpDesignY[key]
+ NewExpDesignY[key] = np.vstack(
+ (oldExpDesignY[key], Y_PC_can[key])
+ )
- MetaModel.ExpDesign.sampling_method = 'user'
- MetaModel.ExpDesign.X = NewExpDesignX
- MetaModel.ExpDesign.Y = NewExpDesignY
+ PCE_Model_can.ExpDesign.sampling_method = 'user'
+ PCE_Model_can.ExpDesign.X = NewExpDesignX
+ PCE_Model_can.ModelOutputDict = NewExpDesignY
+ PCE_Model_can.ExpDesign.Y = NewExpDesignY
# Train the model for the observed data using x_can
- MetaModel.input_obj.poly_coeffs_flag = False
- PCE_Model_can = MetaModel.train_norm_design(Model, parallel=False)
+ PCE_Model_can.input_obj.poly_coeffs_flag = False
+ PCE_Model_can.train_norm_design(parallel=False)
+
+ # Set the ExpDesign to its original values
+ PCE_Model_can.ExpDesign.X = oldExpDesignX
+ PCE_Model_can.ModelOutputDict = oldExpDesignY
+ PCE_Model_can.ExpDesign.Y = oldExpDesignY
if var.lower() == 'mi':
# Mutual information based on Krause et al
@@ -1078,10 +1018,16 @@ class MetaModelEngine():
elif method.lower() == 'bayesactdesign':
NCandidate = candidates.shape[0]
U_J_d = np.zeros((NCandidate))
+ # Evaluate all candidates
+ y_can, std_can = self.MetaModel.eval_metamodel(samples=candidates)
+ # loop through candidates
for idx, X_can in tqdm(enumerate(candidates), ascii=True,
- desc="OptBayesianDesign"):
- U_J_d[idx] = self.util_BayesianActiveDesign(X_can, sigma2Dict,
- var)
+ desc="BAL Design"):
+ y_hat = {key: items[idx] for key, items in y_can.items()}
+ std = {key: items[idx] for key, items in std_can.items()}
+ U_J_d[idx] = self.util_BayesianActiveDesign(
+ y_hat, std, sigma2Dict, var)
+
elif method.lower() == 'bayesoptdesign':
NCandidate = candidates.shape[0]
U_J_d = np.zeros((NCandidate))
@@ -1399,7 +1345,7 @@ class MetaModelEngine():
candidates, n_cand_groups, axis=0
)
- results = Parallel(n_jobs=-1, backend='threading')(
+ results = Parallel(n_jobs=-1, backend='multiprocessing')(
delayed(self.run_util_func)(
exploit_method, split_cand[i], i, sigma2, var, X_MC)
for i in range(n_cand_groups))
@@ -1427,20 +1373,19 @@ class MetaModelEngine():
U_J_d = np.mean(U_J_d.reshape(-1, n_candidates), axis=1)
# create surrogate model for U_J_d
- from sklearn.preprocessing import MinMaxScaler
- # Take care of inf entries
- good_indices = [i for i, arr in enumerate(U_J_d)
- if np.isfinite(arr).all()]
- scaler = MinMaxScaler()
- X_S = scaler.fit_transform(candidates[good_indices])
- gp = MetaModel.gaussian_process_emulator(
- X_S, U_J_d[good_indices], autoSelect=True
- )
- U_J_d = gp.predict(scaler.transform(allCandidates))
+ # from sklearn.preprocessing import MinMaxScaler
+ # # Take care of inf entries
+ # good_indices = [i for i, arr in enumerate(U_J_d)
+ # if np.isfinite(arr).all()]
+ # scaler = MinMaxScaler()
+ # X_S = scaler.fit_transform(candidates[good_indices])
+ # gp = MetaModel.gaussian_process_emulator(
+ # X_S, U_J_d[good_indices], autoSelect=False
+ # )
+ # U_J_d = gp.predict(scaler.transform(allCandidates))
# Normalize U_J_d
norm_U_J_d = U_J_d / np.sum(U_J_d)
- print("norm_U_J_d:\n", norm_U_J_d)
else:
norm_U_J_d = np.zeros((len(scoreExploration)))
@@ -1742,10 +1687,13 @@ class MetaModelEngine():
# Surrogate error if valid dataset is given.
if rmse is not None:
tot_sigma2s += rmse[out]**2
+ else:
+ tot_sigma2s += np.mean(std_pce[out])**2
likelihoods *= stats.multivariate_normal.pdf(
y_hat_pce[out], data, np.diag(tot_sigma2s),
allow_singular=True)
+
self.Likelihoods = likelihoods
return likelihoods