diff --git a/BayesValidRox/surrogate_models/sequential_design.py b/BayesValidRox/surrogate_models/sequential_design.py
new file mode 100644
index 0000000000000000000000000000000000000000..33580f73c30d9272b241dd328d82297993f74dfc
--- /dev/null
+++ b/BayesValidRox/surrogate_models/sequential_design.py
@@ -0,0 +1,1650 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 28 09:21:18 2022
+
+@author: farid
+"""
+import numpy as np
+from scipy import stats, signal, linalg, sparse
+from scipy.spatial import distance
+from copy import deepcopy, copy
+import tqdm
+import scipy.optimize as opt
+from sklearn.metrics import mean_squared_error
+import multiprocessing
+import matplotlib.pyplot as plt
+import sys
+import os
+import seaborn as sns
+
+from .Exploration import Exploration
+
+
+class SeqDesign():
+ def __init__(self, MetaModel):
+ self.MetaModel = MetaModel
+ self.Model = MetaModel.ModelObj
+
+ # -------------------------------------------------------------------------
+ def util_VarBasedDesign(self, X_can, index, UtilMethod='Entropy'):
+ """
+ Computes the exploitation scores based on:
+ active learning MacKay(ALM) and active learning Cohn (ALC)
+ Paper: Sequential Design with Mutual Information for Computer
+ Experiments (MICE): Emulation of a Tsunami Model by Beck and Guillas
+ (2016)
+ """
+ OutputDictY = self.MetaModel.ExpDesign.Y
+ OutputNames = list(OutputDictY.keys())[1:]
+
+ if UtilMethod == 'Entropy':
+ # ----- Entropy/MMSE/active learning MacKay(ALM) -----
+ # Compute perdiction variance of the old model
+ Y_PC_can, std_PC_can = self.MetaModel.eval_metamodel(samples=X_can)
+ canPredVar = {key: std_PC_can[key]**2 for key in OutputNames}
+
+ varPCE = np.zeros((len(OutputNames), X_can.shape[0]))
+ for KeyIdx, key in enumerate(OutputNames):
+ varPCE[KeyIdx] = np.max(canPredVar[key], axis=1)
+ score = np.max(varPCE, axis=0)
+
+ elif UtilMethod == 'EIGF':
+ # ----- Expected Improvement for Global fit -----
+ # Eq (5) from Liu et al.(2018)
+ # Compute perdiction error and variance of the old model
+ 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}
+
+ EIGF_PCE = np.zeros((len(OutputNames), X_can.shape[0]))
+ for KeyIdx, key in enumerate(OutputNames):
+ residual = predError[key] - OutputDictY[key][index]
+ var = canPredVar[key]
+ EIGF_PCE[KeyIdx] = np.max(residual**2 + var, axis=1)
+ score = np.max(EIGF_PCE, axis=0)
+
+ return -1 * score # -1 is for minimization instead of maximization
+
+ # -------------------------------------------------------------------------
+ def util_BayesianActiveDesign(self, X_can, sigma2Dict, var='DKL'):
+
+ # Evaluate the PCE metamodels at that location ???
+ Y_mean_can, Y_std_can = self.eval_metamodel(samples=np.array([X_can]))
+
+ # Get the data
+ obs_data = self.Observations
+ nObs = self.Model.nObs
+ # TODO: Analytical DKL
+ # 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, nObs))
+ logsamples = np.zeros((MCsize, nObs))
+ for i in range(nObs):
+ 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
+
+ # Rejection Step
+ # Random numbers between 0 and 1
+ unif = np.random.rand(1, MCsize)[0]
+
+ # Reject the poorly performed prior
+ accepted = (likelihoods/np.max(likelihoods)) >= unif
+
+ # Prior-based estimation of BME
+ logBME = np.log(np.nanmean(likelihoods))
+
+ # Posterior-based expectation of likelihoods
+ postLikelihoods = likelihoods[accepted]
+ postExpLikelihoods = np.mean(np.log(postLikelihoods))
+
+ # Posterior-based expectation of prior densities
+ postExpPrior = np.mean(logPriorLikelihoods[accepted])
+
+ # Utility function Eq.2 in Ref. (2)
+ # Posterior covariance matrix after observing data y
+ # Kullback-Leibler Divergence (Sergey's paper)
+ if var == 'DKL':
+
+ # TODO: Calculate the correction factor for BME
+ # BMECorrFactor = self.BME_Corr_Weight(PCE_SparseBayes_can,
+ # ObservationData, sigma2Dict)
+ # BME += BMECorrFactor
+ # Haun et al implementation
+ # U_J_d = np.mean(np.log(Likelihoods[Likelihoods!=0])- logBME)
+ U_J_d = postExpLikelihoods - logBME
+
+ # Marginal log likelihood
+ elif var == 'BME':
+ U_J_d = logBME
+
+ # Entropy-based information gain
+ elif var == 'infEntropy':
+ logBME = np.log(np.nanmean(likelihoods))
+ infEntropy = logBME - postExpPrior - postExpLikelihoods
+ U_J_d = infEntropy * -1 # -1 for minimization
+
+ # Bayesian information criterion
+ elif var == 'BIC':
+ coeffs = self.MetaModel.CoeffsDict.values()
+ nModelParams = max(len(v) for val in coeffs for v in val.values())
+ maxL = np.nanmax(likelihoods)
+ U_J_d = -2 * np.log(maxL) + np.log(nObs) * nModelParams
+
+ # Akaike information criterion
+ elif var == 'AIC':
+ coeffs = self.MetaModel.CoeffsDict.values()
+ nModelParams = max(len(v) for val in coeffs for v in val.values())
+ maxlogL = np.log(np.nanmax(likelihoods))
+ AIC = -2 * maxlogL + 2 * nModelParams
+ # 2 * nModelParams * (nModelParams+1) / (nObs-nModelParams-1)
+ penTerm = 0
+ U_J_d = 1*(AIC + penTerm)
+
+ # Deviance information criterion
+ elif var == 'DIC':
+ # 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)
+ DIC = -2 * np.log(Likelihoods_theta_mean) + 2 * N_star_p
+
+ U_J_d = DIC
+
+ else:
+ print('The algorithm you requested has not been implemented yet!')
+
+ # Handle inf and NaN (replace by zero)
+ if np.isnan(U_J_d) or U_J_d == -np.inf or U_J_d == np.inf:
+ U_J_d = 0.0
+
+ return -1 * U_J_d # -1 is for minimization instead of maximization
+
+ # -------------------------------------------------------------------------
+ def util_BayesianDesign(self, X_can, X_MC, sigma2Dict, var='DKL'):
+
+ # To avoid changes ub original aPCE object
+ Model = self.Model
+ PCEModel = deepcopy(self.MetaModel)
+
+ # Old Experimental design
+ oldExpDesignX = PCEModel.ExpDesign.X
+ oldExpDesignY = PCEModel.ExpDesign.Y
+
+ # Evaluate the PCE metamodels at that location ???
+ Y_PC_can, _ = self.MetaModel.eval_metamodel(samples=np.array([X_can]))
+
+ # Add all suggestion as new 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]
+
+ PCEModel.ExpDesign.SamplingMethod = 'user'
+ PCEModel.ExpDesign.X = NewExpDesignX
+ PCEModel.ExpDesign.Y = NewExpDesignY
+
+ # Create the SparseBayes-based PCE metamodel:
+ PCEModel.Inputs.polycoeffsFlag = False
+ univ_p_val = self.MetaModel.univ_basis_vals(X_can)
+ G_n_m_all = np.zeros((len(Model.Output.Names), Model.nObs))
+
+ for idx, key in enumerate(Model.Output.Names):
+ for i in range(Model.nObs):
+ BasisIndices = PCEModel.BasisDict[key]["y_"+str(i+1)]
+ clf_poly = PCEModel.clf_poly[key]["y_"+str(i+1)]
+ Mn = clf_poly.coef_
+ Sn = clf_poly.sigma_
+ beta = clf_poly.alpha_
+ active = clf_poly.active_
+ Psi = self.MetaModel.PCE_create_Psi(BasisIndices, univ_p_val)
+
+ Sn_new_inv = np.linalg.inv(Sn)
+ Sn_new_inv += beta * np.dot(Psi[:, active].T, Psi[:, active])
+ Sn_new = np.linalg.inv(Sn_new_inv)
+
+ Mn_new = np.dot(Sn_new_inv, Mn[active]).reshape(-1, 1)
+ Mn_new += beta * np.dot(Psi[:, active].T, Y_PC_can[key][0, i])
+ Mn_new = np.dot(Sn_new, Mn_new).flatten()
+
+ # Compute new moments
+ mean_old = Mn[0]
+ mean_new = Mn_new[0]
+ std_old = np.sqrt(np.sum(np.square(Mn[1:])))
+ std_new = np.sqrt(np.sum(np.square(Mn_new[1:])))
+
+ G_n_m = np.log(std_old/std_new) - 1./2
+ G_n_m += std_new**2 / (2*std_new**2)
+ G_n_m += (mean_new - mean_old)**2 / (2*std_old**2)
+
+ G_n_m_all[idx, i] = G_n_m
+
+ clf_poly.coef_[active] = Mn_new
+ clf_poly.sigma_ = Sn_new
+ PCEModel.clf_poly[key]["y_"+str(i+1)] = clf_poly
+
+ # return np.sum(G_n_m_all)
+ # PCEModel.create_PCE_normDesign(Model, OptDesignFlag=True)
+ PCE_SparseBayes_can = PCEModel
+
+ # Get the data
+ obs_data = self.MetaModel.Observations
+
+ # ---------- Inner MC simulation for computing Utility Value ----------
+ # Estimation of the integral via Monte Varlo integration
+ MCsize = X_MC.shape[0]
+ ESS = 0
+
+ while ((ESS > MCsize) or (ESS < 1)):
+
+ # Enriching Monte Carlo samples if need be
+ if ESS != 0:
+ X_MC = self.MetaModel.ExpDesign.GetSample(MCsize, 'random')
+
+ # Evaluate the PCEModel at the given samples
+ Y_MC, std_MC = PCE_SparseBayes_can.eval_metamodel(samples=X_MC)
+
+ # 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.sum(likelihoods)))
+
+ # Enlarge sample size if it doesn't fulfill the criteria
+ if ((ESS > MCsize) or (ESS < 1)):
+ MCsize *= 10
+ ESS = 0
+
+ # Rejection Step
+ # Random numbers between 0 and 1
+ unif = np.random.rand(1, MCsize)[0]
+
+ # Reject the poorly performed prior
+ accepted = (likelihoods/np.max(likelihoods)) >= unif
+
+ # -------------------- Utility functions --------------------
+ # Utility function Eq.2 in Ref. (2)
+ # Kullback-Leibler Divergence (Sergey's paper)
+ if var == 'DKL':
+
+ # Prior-based estimation of BME
+ logBME = np.log(np.nanmean(likelihoods))
+
+ # Posterior-based expectation of likelihoods
+ postLikelihoods = likelihoods[accepted]
+ postExpLikelihoods = np.mean(np.log(postLikelihoods))
+
+ # Haun et al implementation
+ U_J_d = np.mean(np.log(likelihoods[likelihoods != 0])- logBME)
+
+ U_J_d = np.sum(G_n_m_all)
+ # Ryan et al (2014) implementation
+ # importanceWeights = Likelihoods[Likelihoods!=0]/np.sum(Likelihoods[Likelihoods!=0])
+ # U_J_d = np.mean(importanceWeights*np.log(Likelihoods[Likelihoods!=0])) - logBME
+
+ # U_J_d = postExpLikelihoods - logBME
+
+ # Marginal likelihood
+ elif var == 'BME':
+
+ # Prior-based estimation of BME
+ logBME = np.log(np.nanmean(likelihoods))
+ U_J_d = logBME
+
+ # Bayes risk likelihood
+ elif var == 'BayesRisk':
+
+ U_J_d = -1 * np.var(likelihoods)
+
+ # Entropy-based information gain
+ elif var == 'infEntropy':
+ # Prior-based estimation of BME
+ logBME = np.log(np.nanmean(likelihoods))
+
+ # Posterior-based expectation of likelihoods
+ postLikelihoods = likelihoods[accepted] / np.nansum(likelihoods[accepted])
+ postExpLikelihoods = np.mean(np.log(postLikelihoods))
+
+ # Posterior-based expectation of prior densities
+ postExpPrior = np.mean(logPriorLikelihoods[accepted])
+
+ infEntropy = logBME - postExpPrior - postExpLikelihoods
+
+ U_J_d = infEntropy * -1 # -1 for minimization
+
+ # D-Posterior-precision
+ elif var == 'DPP':
+ X_Posterior = X_MC[accepted]
+ # covariance of the posterior parameters
+ U_J_d = -np.log(np.linalg.det(np.cov(X_Posterior)))
+
+ # A-Posterior-precision
+ elif var == 'APP':
+ X_Posterior = X_MC[accepted]
+ # trace of the posterior parameters
+ U_J_d = -np.log(np.trace(np.cov(X_Posterior)))
+
+ else:
+ print('The algorithm you requested has not been implemented yet!')
+
+ return -1 * U_J_d # -1 is for minimization instead of maximization
+
+ # -------------------------------------------------------------------------
+ def subdomain(self, Bounds, NrofNewSample):
+
+ NofPa = self.MetaModel.NofPa
+ NofSubdomain = NrofNewSample + 1
+ LinSpace = np.zeros((NofPa, NofSubdomain))
+
+ for i in range(NofPa):
+ LinSpace[i] = np.linspace(start=Bounds[i][0], stop=Bounds[i][1],
+ num=NofSubdomain)
+ Subdomains = []
+ for k in range(NofSubdomain-1):
+ mylist = []
+ for i in range(NofPa):
+ mylist.append((LinSpace[i, k+0], LinSpace[i, k+1]))
+ Subdomains.append(tuple(mylist))
+
+ return Subdomains
+
+ # -------------------------------------------------------------------------
+ def Utility_runner(self, method, allCandidates, index, sigma2Dict=None,
+ var=None, X_MC=None):
+
+ if method == 'VarOptDesign':
+ U_J_d = self.util_VarBasedDesign(allCandidates, index, var)
+
+ elif method == 'BayesActDesign':
+ NCandidate = allCandidates.shape[0]
+ U_J_d = np.zeros((NCandidate))
+ for idx, X_can in tqdm(enumerate(allCandidates), ascii=True,
+ desc="OptBayesianDesign"):
+ U_J_d[idx] = self.util_BayesianActiveDesign(X_can, sigma2Dict,
+ var)
+ elif method == 'BayesOptDesign':
+ NCandidate = allCandidates.shape[0]
+ U_J_d = np.zeros((NCandidate))
+ for idx, X_can in tqdm(enumerate(allCandidates), ascii=True,
+ desc="OptBayesianDesign"):
+ U_J_d[idx] = self.util_BayesianDesign(X_can, X_MC, sigma2Dict,
+ var)
+ return (index, -1 * U_J_d)
+
+ # -------------------------------------------------------------------------
+ def dual_annealing(self, method, Bounds, sigma2Dict, var, Run_No,
+ debugFlag=False):
+
+ Model = self.Model
+ MaxFunItr = self.MetaModel.ExpDesign.MaxFunItr
+
+ if method == 'VarOptDesign':
+ Res_Global = opt.dual_annealing(self.util_VarBasedDesign,
+ bounds=Bounds,
+ args=(Model, var),
+ maxfun=MaxFunItr)
+
+ elif method == 'BayesOptDesign':
+ Res_Global = opt.dual_annealing(self.util_BayesianDesign,
+ bounds=Bounds,
+ args=(Model, sigma2Dict, var),
+ maxfun=MaxFunItr)
+
+ if debugFlag:
+ print(f"global minimum: xmin = {Res_Global.x}, "
+ f"f(xmin) = {Res_Global.fun:.6f}, nfev = {Res_Global.nfev}")
+
+ return (Run_No, Res_Global.x)
+
+ # -------------------------------------------------------------------------
+ def tradoffWeights(self, TradeOffScheme, OldExpDesign, OutputDictY):
+
+ if TradeOffScheme == 'None':
+ weightExploration = 0
+
+ elif TradeOffScheme == 'equal':
+ weightExploration = 0.5
+
+ elif TradeOffScheme == 'epsilon-decreasing':
+ # epsilon-decreasing scheme
+ # Start with more exploration and increase the influence of
+ # exploitation along the way with a exponential decay function
+ initNSamples = self.MetaModel.ExpDesign.initNrSamples
+ maxNSamples = self.MetaModel.ExpDesign.MaxNSamples
+
+ itrNumber = (self.MetaModel.ExpDesign.X.shape[0] - initNSamples)
+ itrNumber //= self.MetaModel.ExpDesign.NrofNewSample
+
+ tau2 = -(maxNSamples-initNSamples-1) / np.log(1e-8)
+ weightExploration = signal.exponential(maxNSamples-initNSamples, 0,
+ tau2, False)[itrNumber]
+
+ elif TradeOffScheme == 'adaptive':
+
+ # Extract itrNumber
+ initNSamples = self.MetaModel.ExpDesign.initNrSamples
+ maxNSamples = self.MetaModel.ExpDesign.MaxNSamples
+ itrNumber = (self.ExpDesign.X.shape[0] - initNSamples)
+ itrNumber //= self.ExpDesign.NrofNewSample
+
+ if itrNumber == 0:
+ weightExploration = 0.5
+ else:
+ # # Extract the model errors from the last and next to last
+ # # iterations
+ # errorModel_i , errorModel_i_1 = self.errorModel[itrNumber],
+ # self.errorModel[itrNumber-1]
+
+ # # Evaluate the error models for all selected samples so far
+ # eLCAllCands_i, _ = errorModel_i.eval_errormodel(OldExpDesign)
+ # eLCAllCands_i_1, _ = errorModel_i_1.eval_errormodel(OldExpDesign)
+
+ # # Local improvement of LC error at last selected design
+ # sl_i = np.max(np.dstack(eLCAllCands_i.values())[-1])
+ # sl_i_1 = np.max(np.dstack(eLCAllCands_i_1.values())[-1])
+
+ # p = sl_i**2 / sl_i_1**2
+
+ # # Global improvement of LC error at OldExpDesign
+ # sg_i = np.max(np.dstack(eLCAllCands_i.values()),axis=1)
+ # sg_i_1 = np.max(np.dstack(eLCAllCands_i_1.values()),axis=1)
+
+ # q = np.sum(np.square(sg_i)) / np.sum(np.square(sg_i_1))
+
+ # weightExploration = min([0.5*p/q, 1])
+
+ # TODO: New adaptive trade-off according to Liu et al. (2017)
+ # Mean squared error for last design point
+ last_EDX = OldExpDesign[-1].reshape(1, -1)
+ lastPCEY, _ = self.MetaModel.eval_metamodel(samples=last_EDX)
+ pce_y = np.array(list(lastPCEY.values()))[:, 0]
+ y = np.array(list(OutputDictY.values())[1:])[:, -1, :]
+ mseError = mean_squared_error(pce_y, y)
+
+ # Mean squared CV - error for last design point
+ error = []
+ for V in self.MetaModel.LCerror.values():
+ for v in V.values():
+ error.append(v[-1])
+ mseCVError = np.mean(np.square(error))
+ weightExploration = 0.99 * min([0.5*mseError/mseCVError, 1])
+
+ return weightExploration
+
+ # -------------------------------------------------------------------------
+ def opt_SeqDesign(self, sigma2Dict, NCandidate=5, var='DKL'):
+
+ # Initialization
+ PCEModel = self.MetaModel
+ Model = self.Model
+ Bounds = PCEModel.BoundTuples
+ NrNewSample = PCEModel.ExpDesign.NrofNewSample
+ ExploreMethod = PCEModel.ExpDesign.ExploreMethod
+ ExploitMethod = PCEModel.ExpDesign.ExploitMethod
+ NrofCandGroups = PCEModel.ExpDesign.NrofCandGroups
+ TradeOffScheme = PCEModel.ExpDesign.TradeOffScheme
+
+ OldExpDesign = PCEModel.ExpDesign.X
+ OutputDictY = PCEModel.ExpDesign.Y.copy()
+ ndim = PCEModel.ExpDesign.X.shape[1]
+
+ # -----------------------------------------
+ # ----------- CUSTOMIZED METHODS ----------
+ # -----------------------------------------
+ # Utility function ExploitMethod provided by user
+ if ExploitMethod.lower() == 'user':
+
+ Xnew, filteredSamples = PCEModel.ExpDesign.ExploitFunction(self)
+
+ print("\n")
+ print("\nXnew:\n", Xnew)
+
+ return Xnew, filteredSamples
+
+ # -----------------------------------------
+ # ---------- EXPLORATION METHODS ----------
+ # -----------------------------------------
+ if ExploreMethod == 'dual annealing':
+ # ------- EXPLORATION: OPTIMIZATION -------
+ import time
+ start_time = time.time()
+
+ # Divide the domain to subdomains
+ args = []
+ Subdomains = self.subdomain(Bounds, NrNewSample)
+ for i in range(NrNewSample):
+ args.append((ExploitMethod, Subdomains[i], sigma2Dict, var, i))
+
+ # Multiprocessing
+ pool = multiprocessing.Pool(multiprocessing.cpu_count())
+
+ # With Pool.starmap_async()
+ results = pool.starmap_async(self.dual_annealing, args).get()
+
+ # Close the pool
+ pool.close()
+
+ Xnew = np.array([results[NofE][1] for NofE in range(NrNewSample)])
+
+ print("\nXnew:\n", Xnew)
+
+ elapsed_time = time.time() - start_time
+ print("\n")
+ print(f"elapsed_time: {round(elapsed_time,2)} sec.")
+ print('-'*20)
+
+ elif ExploreMethod == 'LOOCV':
+ # -----------------------------------------------------------------
+ # TODO: LOOCV model construnction based on Feng et al. (2020)
+ # 'LOOCV':
+ # Initilize the ExploitScore array
+ OutputNames = list(OutputDictY.keys())[1:]
+
+ # Generate random samples
+ allCandidates = PCEModel.ExpDesign.GetSample(
+ NCandidate, SamplingMethod='random'
+ )
+
+ # Construct error model based on LCerror
+ errorModel = PCEModel.create_ModelError(OldExpDesign, self.LCerror)
+ self.errorModel.append(copy(errorModel))
+
+ # Evaluate the error models for allCandidates
+ eLCAllCands, _ = errorModel.eval_errormodel(allCandidates)
+ # Select the maximum as the representative error
+ eLCAllCands = np.dstack(eLCAllCands.values())
+ eLCAllCandidates = np.max(eLCAllCands, axis=1)[:, 0]
+
+ # Normalize the error w.r.t the maximum error
+ scoreExploration = eLCAllCandidates / np.sum(eLCAllCandidates)
+
+ else:
+ # ------- EXPLORATION: SPACE-FILLING DESIGN -------
+ # Generate candidate samples from Exploration class
+ explore = Exploration(PCEModel, NCandidate)
+ explore.w = 100 # * ndim #500
+ # Select criterion (mc-intersite-proj-th, mc-intersite-proj)
+ explore.mcCriterion = 'mc-intersite-proj'
+ allCandidates, scoreExploration = explore.getExplorationSamples()
+
+ # Temp: ---- Plot all candidates -----
+ if ndim == 2:
+ def plotter(points, allCandidates, Method,
+ scoreExploration=None):
+ if Method == 'Voronoi':
+ from scipy.spatial import Voronoi, voronoi_plot_2d
+ vor = Voronoi(points)
+ fig = voronoi_plot_2d(vor)
+ ax1 = fig.axes[0]
+ else:
+ fig = plt.figure()
+ ax1 = fig.add_subplot(111)
+ ax1.scatter(points[:, 0], points[:, 1], s=10, c='r',
+ marker="s", label='Old Design Points')
+ ax1.scatter(allCandidates[:, 0], allCandidates[:, 1], s=10,
+ c='b', marker="o", label='Design candidates')
+ for i in range(points.shape[0]):
+ txt = 'p'+str(i+1)
+ ax1.annotate(txt, (points[i, 0], points[i, 1]))
+ if scoreExploration is not None:
+ for i in range(allCandidates.shape[0]):
+ txt = str(round(scoreExploration[i], 5))
+ ax1.annotate(txt, (allCandidates[i, 0],
+ allCandidates[i, 1]))
+
+ plt.xlim(self.BoundTuples[0])
+ plt.ylim(self.BoundTuples[1])
+ # plt.show()
+ plt.legend(loc='upper left')
+
+ # -----------------------------------------
+ # --------- EXPLOITATION METHODS ----------
+ # -----------------------------------------
+ if ExploitMethod == 'BayesOptDesign' or\
+ ExploitMethod == 'BayesActDesign':
+
+ # ------- Calculate Exoploration weight -------
+ # Compute exploration weight based on trade off scheme
+ weightExploration = self.tradoffWeights(TradeOffScheme,
+ OldExpDesign,
+ OutputDictY)
+ print(f"\nweightExploration={weightExploration:0.3f} "
+ f"weightExploitation={1- weightExploration:0.3f}")
+
+ # ------- EXPLOITATION: BayesOptDesign & ActiveLearning -------
+ if weightExploration != 1.0:
+
+ # Create a sample pool for rejection sampling
+ MCsize = 15000
+ X_MC = PCEModel.ExpDesign.GetSample(MCsize, 'random')
+
+ # Multiprocessing
+ pool = multiprocessing.Pool(multiprocessing.cpu_count())
+
+ # Split the candidates in groups for multiprocessing
+ split_cand = np.array_split(allCandidates,
+ NrofCandGroups, axis=0)
+ args = []
+ for i in range(NrofCandGroups):
+ args.append((ExploitMethod, split_cand[i], i, sigma2Dict,
+ var, X_MC))
+
+ # With Pool.starmap_async()
+ results = pool.starmap_async(self.Utility_runner, args).get()
+
+ # Close the pool
+ pool.close()
+
+ # Retrieve the results and append them
+ U_J_d = np.concatenate([results[NofE][1] for NofE in
+ range(NrofCandGroups)])
+
+ # Get the expected value (mean) of the Utility score
+ # for each cell
+ if ExploreMethod == 'Voronoi':
+ U_J_d = np.mean(U_J_d.reshape(-1, NCandidate), axis=1)
+
+ # 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)))
+
+ # ------- Calculate Total score -------
+ # ------- Trade off between EXPLORATION & EXPLOITATION -------
+ # Total score
+ totalScore = (1 - weightExploration)*norm_U_J_d
+ totalScore += weightExploration*scoreExploration
+
+ # temp: Plot
+ # dim = self.ExpDesign.X.shape[1]
+ # if dim == 2:
+ # plotter(self.ExpDesign.X, allCandidates, ExploreMethod)
+
+ # ------- Select the best candidate -------
+ # find an optimal point subset to add to the initial design by
+ # maximization of the utility score and taking care of NaN values
+ temp = totalScore.copy()
+ temp[np.isnan(totalScore)] = -np.inf
+ sorted_idxtotalScore = np.argsort(temp)[::-1]
+ bestIdx = sorted_idxtotalScore[:NrNewSample]
+
+ # select the requested number of samples
+ if ExploreMethod == 'Voronoi':
+ Xnew = np.zeros((NrNewSample, ndim))
+ for i, idx in enumerate(bestIdx):
+ X_can = explore.closestPoints[idx]
+
+ # Calculate the maxmin score for the region of interest
+ newSamples, maxminScore = explore.getMCSamples(X_can)
+
+ # select the requested number of samples
+ Xnew[i] = newSamples[np.argmax(maxminScore)]
+ else:
+ Xnew = allCandidates[sorted_idxtotalScore[:NrNewSample]]
+
+ elif ExploitMethod == 'VarOptDesign':
+ # ------- EXPLOITATION: VarOptDesign -------
+ UtilMethod = var
+
+ # ------- Calculate Exoploration weight -------
+ # Compute exploration weight based on trade off scheme
+ weightExploration = self.tradoffWeights(TradeOffScheme,
+ OldExpDesign,
+ OutputDictY)
+ print(f"\nweightExploration={weightExploration:0.3f} "
+ f"weightExploitation={1- weightExploration:0.3f}")
+
+ # Generate candidate samples from Exploration class
+ OutputNames = list(OutputDictY.keys())[1:]
+ nMeasurement = OutputDictY[OutputNames[0]].shape[1]
+
+ # Find indices of the Vornoi cells with samples
+ goodSampleIdx = []
+ for idx in range(len(explore.closestPoints)):
+ if len(explore.closestPoints[idx]) != 0:
+ goodSampleIdx.append(idx)
+
+ # Find sensitive region
+ if UtilMethod == 'LOOCV':
+ LCerror = PCEModel.LCerror
+ allModifiedLOO = np.zeros((len(OldExpDesign), len(OutputNames),
+ nMeasurement))
+ for y_idx, y_key in enumerate(OutputNames):
+ for idx, key in enumerate(LCerror[y_key].keys()):
+ allModifiedLOO[:, y_idx, idx] = abs(
+ LCerror[y_key][key])
+
+ ExploitScore = np.max(np.max(allModifiedLOO, axis=1), axis=1)
+
+ elif UtilMethod in ['EIGF', 'Entropy']:
+ # ----- All other in ['EIGF', 'Entropy', 'ALM'] -----
+ # Initilize the ExploitScore array
+ ExploitScore = np.zeros((len(OldExpDesign), len(OutputNames)))
+
+ # Split the candidates in groups for multiprocessing
+ if ExploreMethod != 'Voronoi':
+ split_cand = np.array_split(allCandidates,
+ NrofCandGroups,
+ axis=0)
+ goodSampleIdx = range(NrofCandGroups)
+ else:
+ split_cand = explore.closestPoints
+
+ # Split the candidates in groups for multiprocessing
+ args = []
+ for index in goodSampleIdx:
+ args.append((ExploitMethod, split_cand[index], index,
+ sigma2Dict, var))
+
+ # Multiprocessing
+ pool = multiprocessing.Pool(multiprocessing.cpu_count())
+ # With Pool.starmap_async()
+ results = pool.starmap_async(self.Utility_runner, args).get()
+
+ # Close the pool
+ pool.close()
+
+ # Retrieve the results and append them
+ ExploitScore = np.concatenate([results[NofE][1] for NofE in
+ range(len(goodSampleIdx))])
+
+ else:
+ raise NameError('The requested utility function is not '
+ 'available.')
+
+ # print("ExploitScore:\n", ExploitScore)
+
+ # find an optimal point subset to add to the initial design by
+ # maximization of the utility score and taking care of NaN values
+ # Total score
+ # Normalize U_J_d
+ ExploitScore = ExploitScore / np.sum(ExploitScore)
+ totalScore = (1 - weightExploration) * ExploitScore
+ totalScore += weightExploration * scoreExploration
+
+ temp = totalScore.copy()
+ sorted_idxtotalScore = np.argsort(temp, axis=0)[::-1]
+ bestIdx = sorted_idxtotalScore[:NrNewSample]
+
+ Xnew = np.zeros((NrNewSample, ndim))
+ if ExploreMethod != 'Voronoi':
+ Xnew = allCandidates[bestIdx]
+ else:
+ for i, idx in enumerate(bestIdx.flatten()):
+ X_can = explore.closestPoints[idx]
+ # plotter(self.ExpDesign.X, X_can, ExploreMethod,
+ # scoreExploration=None)
+
+ # Calculate the maxmin score for the region of interest
+ newSamples, maxminScore = explore.getMCSamples(X_can)
+
+ # select the requested number of samples
+ Xnew[i] = newSamples[np.argmax(maxminScore)]
+
+ elif ExploitMethod == 'alphabetic':
+ # ------- EXPLOITATION: ALPHABETIC -------
+ Xnew = self.util_AlphOptDesign(allCandidates, var)
+
+ elif ExploitMethod == 'Space-filling':
+ # ------- EXPLOITATION: SPACE-FILLING -------
+ totalScore = scoreExploration
+
+ # ------- Select the best candidate -------
+ # find an optimal point subset to add to the initial design by
+ # maximization of the utility score and taking care of NaN values
+ temp = totalScore.copy()
+ temp[np.isnan(totalScore)] = -np.inf
+ sorted_idxtotalScore = np.argsort(temp)[::-1]
+
+ # select the requested number of samples
+ Xnew = allCandidates[sorted_idxtotalScore[:NrNewSample]]
+
+ else:
+ raise NameError('The requested design method is not available.')
+
+ print("\n")
+ print("\nRun No. {}:".format(OldExpDesign.shape[0]+1))
+ print("Xnew:\n", Xnew)
+
+ return Xnew, None
+
+ # -------------------------------------------------------------------------
+ def util_AlphOptDesign(self, allCandidates, var='D-Opt'):
+ """
+ Enrich the Experimental design with the requested alphabetic criterion
+ based on exploring the space with number of sampling points.
+
+ Ref: Hadigol, M., & Doostan, A. (2018). Least squares polynomial chaos
+ expansion: A review of sampling strategies., Computer Methods in
+ Applied Mechanics and Engineering, 332, 382-407.
+
+ Arguments
+ ---------
+ NCandidate : int
+ Number of candidate points to be searched
+
+ var : string
+ Alphabetic optimality criterion
+
+ Returns
+ -------
+ X_new : ndarray[1, n]
+ The new sampling location in the input space.
+ """
+ PCEModelOrig = self.PCEModel
+ Model = self.ModelObj
+ NrofNewSample = PCEModelOrig.ExpDesign.NrofNewSample
+ index = PCEModelOrig.index
+ NCandidate = allCandidates.shape[0]
+
+ # TODO: Loop over outputs
+ OutputName = Model.Output.Names[0]
+
+ # To avoid changes ub original aPCE object
+ PCEModel = deepcopy(PCEModelOrig)
+
+ # Old Experimental design
+ oldExpDesignX = PCEModel.ExpDesign.X
+
+ # TODO: Only one psi can be selected.
+ # Suggestion: Go for the one with the highest LOO error
+ Scores = list(PCEModel.ScoresDict[OutputName].values())
+ ModifiedLOO = [1-score for score in Scores]
+ outIdx = np.argmax(ModifiedLOO[:index])
+
+ # Initialize Phi to save the criterion's values
+ Phi = np.zeros((NCandidate))
+
+ BasisIndices = PCEModelOrig.BasisDict[OutputName]["y_"+str(outIdx+1)]
+ P = len(BasisIndices)
+
+ # ------ Old Psi ------------
+ univ_p_val = PCEModelOrig.univ_basis_vals(oldExpDesignX)
+ Psi = PCEModelOrig.PCE_create_Psi(BasisIndices, univ_p_val)
+
+ # ------ New candidates (Psi_c) ------------
+ # Assemble Psi_c
+ univ_p_val_c = self.univ_basis_vals(allCandidates)
+ Psi_c = self.PCE_create_Psi(BasisIndices, univ_p_val_c)
+
+ for idx in range(NCandidate):
+
+ # Include the new row to the original Psi
+ Psi_cand = np.vstack((Psi, Psi_c[idx]))
+
+ # Information matrix
+ PsiTPsi = np.dot(Psi_cand.T, Psi_cand)
+ M = PsiTPsi / (len(oldExpDesignX)+1)
+
+ if np.linalg.cond(PsiTPsi) > 1e-12 \
+ and np.linalg.cond(PsiTPsi) < 1 / sys.float_info.epsilon:
+ # faster
+ invM = linalg.solve(M, sparse.eye(PsiTPsi.shape[0]).toarray())
+ else:
+ # stabler
+ invM = np.linalg.pinv(M)
+
+ # ---------- Calculate optimality criterion ----------
+ # Optimality criteria according to Section 4.5.1 in Ref.
+
+ # D-Opt
+ if var == 'D-Opt':
+ Phi[idx] = (np.linalg.det(invM)) ** (1/P)
+
+ # A-Opt
+ elif var == 'A-Opt':
+ Phi[idx] = np.trace(invM)
+
+ # K-Opt
+ elif var == 'K-Opt':
+ Phi[idx] = np.linalg.cond(M)
+
+ else:
+ raise Exception('The optimality criterion you requested has '
+ 'not been implemented yet!')
+
+ # find an optimal point subset to add to the initial design
+ # by minimization of the Phi
+ sorted_idxtotalScore = np.argsort(Phi)
+
+ # select the requested number of samples
+ Xnew = allCandidates[sorted_idxtotalScore[:NrofNewSample]]
+
+ return Xnew
+
+ # -------------------------------------------------------------------------
+ def normpdf(self, PCEOutputs, std_PC_MC, obs_data, Sigma2s):
+
+ output_names = self.Model.Output.Names
+
+ SampleSize, index = PCEOutputs[output_names[0]].shape
+ if type(self.MetaModel.index) is list:
+ index = self.MetaModel.index
+ else:
+ index = [self.MetaModel.index]*len(output_names)
+
+ # Flatten the ObservationData
+ TotalData = obs_data[output_names].to_numpy().flatten('F')
+
+ # Remove NaN
+ TotalData = TotalData[~np.isnan(TotalData)]
+ Sigma2s = Sigma2s[~np.isnan(Sigma2s)]
+
+ # Flatten the Output
+ TotalOutputs = np.empty((SampleSize, 0))
+ for idx, key in enumerate(output_names):
+ TotalOutputs = np.hstack((TotalOutputs,
+ PCEOutputs[key][:, :index[idx]]))
+
+ # Covariance Matrix
+ covMatrix = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
+ np.fill_diagonal(covMatrix, Sigma2s)
+
+ # Add the std of the PCE.
+ covMatrix_PCE = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
+ stdPCE = np.empty((SampleSize, 0))
+ for idx, key in enumerate(output_names):
+ stdPCE = np.hstack((stdPCE, std_PC_MC[key][:, :index[idx]]))
+
+ # Expected value of variance (Assump: i.i.d stds)
+ varPCE = np.mean(stdPCE**2, axis=0)
+ # # varPCE = np.var(stdPCE, axis=1)
+ np.fill_diagonal(covMatrix_PCE, varPCE)
+
+ # Aggregate the cov matrices
+ covMatrix += covMatrix_PCE
+
+ # Compute likelihood
+ self.Likelihoods = stats.multivariate_normal.pdf(TotalOutputs,
+ mean=TotalData,
+ cov=covMatrix,
+ allow_singular=True)
+ return self.Likelihoods
+
+ # -------------------------------------------------------------------------
+ def posteriorPlot(self, Posterior, MAP, parNames, key, figsize=(10, 10)):
+
+ # Initialization
+ newpath = (r'Outputs_SeqPosteriorComparison')
+ if not os.path.exists(newpath):
+ os.makedirs(newpath)
+
+ lw = 3.
+ BoundTuples = self.BoundTuples
+ NofPa = len(MAP)
+
+ # This is the true mean of the second mode that we used above:
+ value1 = MAP
+
+ # This is the empirical mean of the sample:
+ value2 = np.mean(Posterior, axis=0)
+
+ if NofPa == 2:
+
+ figPosterior, ax = plt.subplots()
+ plt.hist2d(Posterior[:, 0], Posterior[:, 1], bins=(200, 200),
+ range=np.array([BoundTuples[0], BoundTuples[1]]),
+ cmap=plt.cm.jet)
+
+ plt.xlabel(parNames[0])
+ plt.ylabel(parNames[1])
+
+ plt.xlim(BoundTuples[0])
+ plt.ylim(BoundTuples[1])
+
+ ax.axvline(value1[0], color="g", lw=lw)
+ ax.axhline(value1[1], color="g", lw=lw)
+ ax.plot(value1[0], value1[1], "sg", lw=lw+1)
+
+ ax.axvline(value2[0], ls='--', color="r", lw=lw)
+ ax.axhline(value2[1], ls='--', color="r", lw=lw)
+ ax.plot(value2[0], value2[1], "sr", lw=lw+1)
+
+ else:
+ import corner
+ figPosterior = corner.corner(Posterior, labels=parNames,
+ title_fmt='.2e', show_titles=True,
+ title_kwargs={"fontsize": 12})
+
+ # Extract the axes
+ axes = np.array(figPosterior.axes).reshape((NofPa, NofPa))
+
+ # Loop over the diagonal
+ for i in range(NofPa):
+ ax = axes[i, i]
+ ax.axvline(value1[i], color="g")
+ ax.axvline(value2[i], ls='--', color="r")
+
+ # Loop over the histograms
+ for yi in range(NofPa):
+ for xi in range(yi):
+ ax = axes[yi, xi]
+ ax.axvline(value1[xi], color="g")
+ ax.axvline(value2[xi], ls='--', color="r")
+ ax.axhline(value1[yi], color="g")
+ ax.axhline(value2[yi], ls='--', color="r")
+ ax.plot(value1[xi], value1[yi], "sg")
+ ax.plot(value2[xi], value2[yi], "sr")
+
+ figPosterior.savefig(f'./{newpath}/{key}.svg', bbox_inches='tight')
+ plt.close()
+
+ # Save the posterior as .npy
+ np.save(f'./{newpath}/{key}.npy', Posterior)
+
+ return figPosterior
+
+ # -------------------------------------------------------------------------
+ def hellinger_distance(self, P, Q):
+ """
+ Hellinger distance between two continuous distributions.
+
+ The maximum distance 1 is achieved when P assigns probability zero to
+ every set to which Q assigns a positive probability, and vice versa.
+ 0 (identical) and 1 (maximally different)
+
+ Parameters
+ ----------
+ P : array
+ Reference likelihood.
+ Q : array
+ Estimated likelihood.
+
+ Returns
+ -------
+ float
+ Hellinger distance of two distributions.
+
+ """
+ mu1 = P.mean()
+ Sigma1 = np.std(P)
+
+ mu2 = Q.mean()
+ Sigma2 = np.std(Q)
+
+ term1 = np.sqrt(2*Sigma1*Sigma2 / (Sigma1**2 + Sigma2**2))
+
+ term2 = np.exp(-.25 * (mu1 - mu2)**2 / (Sigma1**2 + Sigma2**2))
+
+ H_squared = 1 - term1 * term2
+
+ return np.sqrt(H_squared)
+
+ # -------------------------------------------------------------------------
+ def BME_Calculator(self, PCEModel, obs_data, sigma2Dict):
+ """
+ This function computes the Bayesian model evidence (BME) via Monte
+ Carlo integration.
+
+ """
+ # Initializations
+ validLikelihoods = PCEModel.validLikelihoods
+ PostSnapshot = PCEModel.ExpDesign.PostSnapshot
+ if PostSnapshot or len(validLikelihoods) != 0:
+ newpath = (r'Outputs_SeqPosteriorComparison')
+ if not os.path.exists(newpath):
+ os.makedirs(newpath)
+
+ SamplingMethod = 'random'
+ MCsize = 100000
+ ESS = 0
+
+ # Estimation of the integral via Monte Varlo integration
+ while ((ESS > MCsize) or (ESS < 1)):
+
+ # Generate samples for Monte Carlo simulation
+ if len(validLikelihoods) == 0:
+ X_MC = PCEModel.ExpDesign.GetSample(MCsize, SamplingMethod)
+ else:
+ X_MC = PCEModel.validSamples
+ MCsize = X_MC.shape[0]
+
+ # Monte Carlo simulation for the candidate design
+ Y_MC, std_MC = PCEModel.eval_metamodel(samples=X_MC)
+
+ # 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 (1000<ESS<MCsize)
+ ESS = 1 / np.sum(np.square(Likelihoods/np.sum(Likelihoods)))
+
+ # Enlarge sample size if it doesn't fulfill the criteria
+ if ((ESS > MCsize) or (ESS < 1)):
+ print('ESS={0} MC size should be larger.'.format(ESS))
+ MCsize = MCsize * 10
+ ESS = 0
+
+ # Rejection Step
+ # Random numbers between 0 and 1
+ unif = np.random.rand(1, MCsize)[0]
+
+ # Reject the poorly performed prior
+ accepted = (Likelihoods/np.max(Likelihoods)) >= unif
+ X_Posterior = X_MC[accepted]
+
+ # ------------------------------------------------------------
+ # --- Kullback-Leibler Divergence & Information Entropy ------
+ # ------------------------------------------------------------
+ # Prior-based estimation of BME
+ logBME = np.log(np.nanmean(Likelihoods))
+
+ # Posterior-based expectation of likelihoods
+ postExpLikelihoods = np.mean(np.log(Likelihoods[accepted]))
+
+ # Posterior-based expectation of prior densities
+ # postExpPrior = np.mean([log_prior(sample) for sample in X_Posterior])
+
+ # Calculate Kullback-Leibler Divergence
+ # KLD = np.mean(np.log(Likelihoods[Likelihoods!=0])- logBME)
+ KLD = postExpLikelihoods - logBME
+
+ # Information Entropy based on Entropy paper Eq. 38
+ # infEntropy = logBME - postExpPrior - postExpLikelihoods
+
+ # If PostSnapshot is True, plot likelihood vs refrence
+ if PostSnapshot or len(validLikelihoods) != 0:
+ idx = len([name for name in os.listdir(newpath) if 'Likelihoods_'
+ in name and os.path.isfile(os.path.join(newpath, name))])
+ fig, ax = plt.subplots()
+ sns.kdeplot(np.log(validLikelihoods[validLikelihoods > 0]),
+ shade=True, color="g", label='Ref. Likelihood')
+ sns.kdeplot(np.log(Likelihoods[Likelihoods > 0]), shade=True,
+ color="b", label='Likelihood with PCE')
+
+ # Hellinger distance
+ ref_like = np.log(validLikelihoods[validLikelihoods > 0])
+ est_like = np.log(Likelihoods[Likelihoods > 0])
+ distHellinger = self.hellinger_distance(ref_like, est_like)
+ text = f"Hellinger Dist.={distHellinger:.3f}\n logBME={logBME:.3f}"
+ "\n DKL={KLD:.3f}"
+
+ plt.text(0.05, 0.75, text, bbox=dict(facecolor='wheat',
+ edgecolor='black',
+ boxstyle='round,pad=1'),
+ transform=ax.transAxes)
+
+ fig.savefig(f'./{newpath}/Likelihoods_{idx}.svg',
+ bbox_inches='tight')
+ plt.close()
+
+ else:
+ distHellinger = 0.0
+
+ return (logBME, KLD, X_Posterior, distHellinger)
+
+ # -------------------------------------------------------------------------
+ def validError_(self):
+
+ PCEModel = self.MetaModel
+ Model = self.Model
+ OutputName = Model.Output.Names
+
+ # Generate random samples
+ Samples = PCEModel.validSamples
+
+ # Extract the original model with the generated samples
+ ModelOutputs = PCEModel.validModelRuns
+
+ # Run the PCE model with the generated samples
+ PCEOutputs, PCEOutputs_std = PCEModel.eval_metamodel(samples=Samples)
+
+ validError_dict = {}
+ # Loop over the keys and compute RMSE error.
+ for key in OutputName:
+ weight = np.mean(np.square(PCEOutputs_std[key]), axis=0)
+ if all(weight == 0):
+ weight = 'variance_weighted'
+ validError_dict[key] = mean_squared_error(ModelOutputs[key],
+ PCEOutputs[key])
+ validError_dict[key] /= np.var(ModelOutputs[key], ddof=1)
+
+ return validError_dict
+
+ # -------------------------------------------------------------------------
+ def error_Mean_Std(self):
+
+ PCEModel = self.MetaModel
+ # Extract the mean and std provided by user
+ df_MCReference = PCEModel.ModelObj.MCReference
+
+ # Compute the mean and std based on the PCEModel
+ PCEMeans = dict()
+ PCEStds = dict()
+ for Outkey, ValuesDict in PCEModel.CoeffsDict.items():
+ PCEMean = np.zeros((len(ValuesDict)))
+ PCEStd = np.zeros((len(ValuesDict)))
+
+ for Inkey, InIdxValues in ValuesDict.items():
+ idx = int(Inkey.split('_')[1]) - 1
+ coeffs = PCEModel.CoeffsDict[Outkey][Inkey]
+
+ # Mean = c_0
+ if coeffs[0] != 0:
+ PCEMean[idx] = coeffs[0]
+ else:
+ PCEMean[idx] = PCEModel.clf_poly[Outkey][Inkey].intercept_
+
+ # Std = sqrt(sum(coeffs[1:]**2))
+ PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
+
+ if PCEModel.DimRedMethod.lower() == 'pca':
+ PCA = PCEModel.pca[Outkey]
+ PCEMean = PCA.mean_ + np.dot(PCEMean, PCA.components_)
+ PCEStd = np.dot(PCEStd, PCA.components_)
+
+ # Compute the error between mean and std of PCEModel and OrigModel
+ RMSE_Mean = mean_squared_error(df_MCReference['mean'], PCEMean,
+ squared=False)
+ RMSE_std = mean_squared_error(df_MCReference['std'], PCEStd,
+ squared=False)
+
+ PCEMeans[Outkey] = PCEMean
+ PCEStds[Outkey] = PCEStd
+
+ return RMSE_Mean, RMSE_std
+
+ # -------------------------------------------------------------------------
+ def create_PCE_SeqDesign(self, Model):
+
+ # Initialization
+ PCEModel = self.MetaModel
+ errorIncreases = False
+ PCEModel.SeqModifiedLOO = {}
+ PCEModel.seqValidError = {}
+ PCEModel.SeqBME = {}
+ PCEModel.SeqKLD = {}
+ PCEModel.SeqDistHellinger = {}
+ PCEModel.seqRMSEMean = {}
+ PCEModel.seqRMSEStd = {}
+ PCEModel.seqMinDist = []
+ pce = True if PCEModel.metaModel.lower() != 'gpe' else False
+ mc_ref = True if hasattr(Model, 'MCReference') else False
+ if mc_ref:
+ Model.read_MCReference()
+
+ # Get the parameters
+ max_n_samples = PCEModel.ExpDesign.MaxNSamples
+ mod_LOO_threshold = PCEModel.ExpDesign.ModifiedLOOThreshold
+ n_canddidate = PCEModel.ExpDesign.NCandidate
+ post_snapshot = PCEModel.ExpDesign.PostSnapshot
+ n_replication = PCEModel.ExpDesign.nReprications
+
+ # Handle if only one UtilityFunctions is provided
+ if not isinstance(PCEModel.ExpDesign.UtilityFunction, list):
+ util_func = [PCEModel.ExpDesign.UtilityFunction]
+
+ # ---------- Initial PCEModel ----------
+ PCEModel = PCEModel.create_PCE_normDesign(Model)
+ initPCEModel = deepcopy(PCEModel)
+
+ # TODO: Loop over outputs
+ OutputName = Model.Output.Names
+
+ # Estimation of the integral via Monte Varlo integration
+ obs_data = PCEModel.Observations
+
+ # Check if data is provided
+ TotalSigma2 = np.empty((0, 1))
+ if len(obs_data) != 0:
+ # ------ Prepare diagonal enteries for co-variance matrix ---------
+ for keyIdx, key in enumerate(Model.Output.Names):
+ # optSigma = 'B'
+ sigma2 = np.array(PCEModel.Discrepancy.Parameters[key])
+ TotalSigma2 = np.append(TotalSigma2, sigma2)
+
+ # Calculate the initial BME
+ out = self.BME_Calculator(initPCEModel, obs_data, TotalSigma2)
+ initBME, initKLD, initPosterior, initDistHellinger = out
+ print("\nInitial BME:", initBME)
+ print("Initial KLD:", initKLD)
+
+ # Posterior snapshot (initial)
+ if post_snapshot:
+ MAP = PCEModel.ExpDesign.MAP
+ parNames = PCEModel.ExpDesign.parNames
+ print('Posterior snapshot (initial) is being plotted...')
+ self.posteriorPlot(initPosterior, MAP, parNames,
+ 'SeqPosterior_init')
+
+ # Check the convergence of the Mean & Std
+ if mc_ref and pce:
+ initRMSEMean, initRMSEStd = self.error_Mean_Std()
+ print(f"Initial Mean and Std error: {initRMSEMean}, {initRMSEStd}")
+
+ # Read the initial experimental design
+ Xinit = initPCEModel.ExpDesign.X
+ initTotalNSamples = len(PCEModel.ExpDesign.X)
+ initYprev = initPCEModel.ModelOutputDict
+ initLCerror = initPCEModel.LCerror
+
+ # Read the initial ModifiedLOO
+ if pce:
+ Scores_all, varExpDesignY = [], []
+ for OutName in OutputName:
+ y = initPCEModel.ExpDesign.Y[OutName]
+ Scores_all.append(list(
+ initPCEModel.ScoresDict[OutName].values()))
+ if PCEModel.DimRedMethod.lower() == 'pca':
+ pca = PCEModel.pca[OutName]
+ components = pca.transform(y)
+ varExpDesignY.append(np.var(components, axis=0))
+ else:
+ varExpDesignY.append(np.var(y, 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)]
+
+ if len(PCEModel.validModelRuns) != 0:
+ initValidError = self.validError_()
+ initValidError = list(initValidError.values())
+ print("\nInitial ValidError:", initValidError)
+
+ # Replicate the sequential design
+ for repIdx in range(n_replication):
+ print(f'>>>> Replication: {repIdx+1}<<<<')
+
+ # To avoid changes ub original aPCE object
+ # PCEModel = copy.deepcopy(initPCEModel)
+ PCEModel.ExpDesign.X = Xinit
+ PCEModel.ExpDesign.Y = initYprev
+ PCEModel.LCerror = initLCerror
+
+ for util_f in util_func:
+ print(f'>>>> UtilityFunction: {util_f} <<<<')
+ # To avoid changes ub original aPCE object
+ # PCEModel = copy.deepcopy(initPCEModel)
+ PCEModel.ExpDesign.X = Xinit
+ PCEModel.ExpDesign.Y = initYprev
+ PCEModel.LCerror = initLCerror
+
+ # Set the experimental design
+ Xprev = Xinit
+ total_n_samples = initTotalNSamples
+ Yprev = initYprev
+
+ Xfull = []
+ Yfull = []
+
+ # Store the initial ModifiedLOO
+ if pce:
+ print("\nInitial ModifiedLOO:", initModifiedLOO)
+ ModifiedLOO = initModifiedLOO
+ SeqModifiedLOO = np.array(ModifiedLOO)
+
+ if len(PCEModel.validModelRuns) != 0:
+ ValidError = initValidError
+ SeqValidError = np.array(ValidError)
+
+ # Check if data is provided
+ if len(obs_data) != 0:
+ SeqBME = np.array([initBME])
+ SeqKLD = np.array([initKLD])
+ SeqDistHellinger = np.array([initDistHellinger])
+
+ if mc_ref and pce:
+ seqRMSEMean = np.array([initRMSEMean])
+ seqRMSEStd = np.array([initRMSEStd])
+
+ # Start Sequential Experimental Design
+ postcnt = 1
+ itrNr = 1
+ while total_n_samples < max_n_samples:
+
+ # Optimal Bayesian Design
+ PCEModel.ExpDesignFlag = 'sequential'
+ Xnew, updatedPrior = self.opt_SeqDesign(TotalSigma2,
+ n_canddidate,
+ util_f)
+
+ S = np.min(distance.cdist(Xinit, Xnew, 'euclidean'))
+ PCEModel.seqMinDist.append(S)
+ print("\nmin Dist from OldExpDesign:", S)
+ print("\n")
+
+ # Evaluate the full model response at the new sample:
+ Ynew, _ = Model.Run_Model_Parallel(
+ Xnew, prevRun_No=total_n_samples
+ )
+ total_n_samples += Xnew.shape[0]
+
+ # ------ Plot the surrogate model vs Origninal Model ------
+ if PCEModel.adaptVerbose:
+ from PostProcessing.adaptPlot import adaptPlot
+ y_hat, std_hat = PCEModel.eval_metamodel(samples=Xnew)
+ adaptPlot(PCEModel, Ynew, y_hat, std_hat, plotED=False)
+
+ # -------- Retrain the surrogate model -------
+ # Extend new experimental design
+ Xfull = np.vstack((Xprev, Xnew))
+
+ # Updating existing key's value
+ for OutIdx in range(len(OutputName)):
+ OutName = OutputName[OutIdx]
+ try:
+ Yfull = np.vstack((Yprev[OutName], Ynew[OutName]))
+ except:
+ Yfull = np.vstack((Yprev[OutName], Ynew))
+
+ PCEModel.ModelOutputDict[OutName] = Yfull
+
+ PCEModel.ExpDesign.SamplingMethod = 'user'
+ PCEModel.ExpDesign.X = Xfull
+ PCEModel.ExpDesign.Y = PCEModel.ModelOutputDict
+
+ # save the Experimental Design for next iteration
+ Xprev = Xfull
+ Yprev = PCEModel.ModelOutputDict
+
+ # Pass the new prior as the input
+ PCEModel.Inputs.polycoeffsFlag = False
+ if updatedPrior is not None:
+ PCEModel.Inputs.polycoeffsFlag = True
+ print("updatedPrior:", updatedPrior.shape)
+ # Arbitrary polynomial chaos
+ for i in range(updatedPrior.shape[1]):
+ PCEModel.Inputs.Marginals[i].DistType = None
+ x = updatedPrior[:, i]
+ PCEModel.Inputs.Marginals[i].InputValues = x
+
+ prevPCEModel = PCEModel
+ PCEModel = PCEModel.create_PCE_normDesign(Model)
+
+ # -------- Evaluate the retrain surrogate model -------
+ # Compute the validation error
+ if len(PCEModel.validModelRuns) != 0:
+ validError = self.validError_()
+ ValidError = list(validError.values())
+ print("\nUpdated ValidError:", ValidError)
+
+ # Extract Modified LOO from Output
+ if pce:
+ Scores_all, varExpDesignY = [], []
+ for OutName in OutputName:
+ y = initPCEModel.ExpDesign.Y[OutName]
+ Scores_all.append(list(
+ PCEModel.ScoresDict[OutName].values()))
+ if PCEModel.DimRedMethod.lower() == 'pca':
+ pca = PCEModel.pca[OutName]
+ components = pca.transform(y)
+ varExpDesignY.append(np.var(components,
+ axis=0))
+ else:
+ varExpDesignY.append(np.var(y, 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], weights=weights)]
+
+ print('\n')
+ print(f"Updated ModifiedLOO {util_f}:\n", ModifiedLOO)
+ print("Xfull:", Xfull.shape)
+ print('\n')
+
+ # check the direction of the error (on average):
+ # if it increases consistently stop the iterations
+ n_checks = 3
+ if itrNr > n_checks * PCEModel.ExpDesign.NrofNewSample:
+ # ss<0 error increasing
+ ss = np.sign(SeqModifiedLOO - ModifiedLOO)
+ errorIncreases = np.sum(np.mean(ss[-2:], axis=1)) <= \
+ -1*n_checks
+
+ # If error is increasing in the last n_check iteration,
+ # stop the search and return the previous PCEModel
+ if errorIncreases:
+ print("Warning: The modified error is increasing "
+ "compared to the last {n_checks} iterations.")
+ PCEModel = prevPCEModel
+ break
+ else:
+ prevPCEModel = PCEModel
+
+ # Store updated ModifiedLOO
+ if pce:
+ SeqModifiedLOO = np.vstack(
+ (SeqModifiedLOO, ModifiedLOO))
+ if len(PCEModel.validModelRuns) != 0:
+ SeqValidError = np.vstack(
+ (SeqValidError, ValidError))
+
+ # -------- Caclulation of BME as accuracy metric -------
+ # Check if data is provided
+ if len(obs_data) != 0:
+ # Calculate the initial BME
+ out = self.BME_Calculator(PCEModel, obs_data,
+ TotalSigma2)
+ BME, KLD, Posterior, DistHellinger = out
+ print('\n')
+ print("Updated BME:", BME)
+ print("Updated KLD:", KLD)
+ print('\n')
+
+ # Plot some snapshots of the posterior
+ stepSnapshot = PCEModel.ExpDesign.stepSnapshot
+ if post_snapshot and postcnt % stepSnapshot == 0:
+ MAP = PCEModel.ExpDesign.MAP
+ parNames = PCEModel.ExpDesign.parNames
+ print('Posterior snapshot is being plotted...')
+ self.posteriorPlot(Posterior, MAP, parNames,
+ 'SeqPosterior_{postcnt}')
+ postcnt += 1
+
+ # Check the convergence of the Mean&Std
+ if mc_ref and pce:
+ print('\n')
+ RMSE_Mean, RMSE_std = self.error_Mean_Std()
+ print(f"Updated Mean and Std error: {RMSE_Mean}, "
+ f"{RMSE_std}")
+ print('\n')
+
+ # Store the updated BME & KLD
+ # Check if data is provided
+ if len(obs_data) != 0:
+ SeqBME = np.vstack((SeqBME, BME))
+ SeqKLD = np.vstack((SeqKLD, KLD))
+ SeqDistHellinger = np.vstack((SeqDistHellinger,
+ DistHellinger))
+ if mc_ref and pce:
+ seqRMSEMean = np.vstack((seqRMSEMean, RMSE_Mean))
+ seqRMSEStd = np.vstack((seqRMSEStd, RMSE_std))
+
+ if pce and any(LOO < mod_LOO_threshold
+ for LOO in ModifiedLOO):
+ break
+ itrNr += 1
+ # Store updated ModifiedLOO and BME in dictonary
+ strKey = f'{util_f}_rep_{repIdx+1}'
+ if pce:
+ PCEModel.SeqModifiedLOO[strKey] = SeqModifiedLOO
+ if len(PCEModel.validModelRuns) != 0:
+ PCEModel.seqValidError[strKey] = SeqValidError
+
+ # Check if data is provided
+ if len(obs_data) != 0:
+ PCEModel.SeqBME[strKey] = SeqBME
+ PCEModel.SeqKLD[strKey] = SeqKLD
+ if len(PCEModel.validLikelihoods) != 0:
+ PCEModel.SeqDistHellinger[strKey] = SeqDistHellinger
+ if mc_ref and pce:
+ PCEModel.seqRMSEMean[strKey] = seqRMSEMean
+ PCEModel.seqRMSEStd[strKey] = seqRMSEStd
+
+ return PCEModel
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 080e09a079d5b16a682515693f1d62a44807bdbd..5511d0afe39680702502dcaa4f0d1f32f4a75c5f 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -5,109 +5,86 @@ Created on Sat Aug 24 13:07:07 2019
@author: farid
"""
-
-
-# To supress warnings
import warnings
-warnings.filterwarnings("ignore")
import numpy as np
import math
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))
-plt.rc('font', family='serif', serif='Arial')
-plt.rc('font', size=SIZE)
-plt.rc('axes', grid = True)
-plt.rc('text', usetex=True)
-plt.rc('axes', linewidth=3)
-plt.rc('axes', grid=True)
-plt.rc('grid', linestyle="-")
-plt.rc('axes', titlesize=SIZE) # fontsize of the axes title
-plt.rc('axes', labelsize=SIZE) # fontsize of the x and y labels
-plt.rc('xtick', labelsize=SIZE) # fontsize of the tick labels
-plt.rc('ytick', labelsize=SIZE) # fontsize of the tick labels
-plt.rc('legend', fontsize=SIZE) # legend fontsize
-plt.rc('figure', titlesize=SIZE) # fontsize of the figure title
-
import scipy as sp
from tqdm import tqdm
-import warnings
-from itertools import combinations, product
-from scipy import stats
+from sklearn.decomposition import PCA as sklearnPCA
import sklearn.linear_model as lm
-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
-import copy
-import seaborn as sns
-import pandas as pd
-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 joblib import Parallel, delayed
+
from .ExpDesigns import ExpDesigns
from .glexindex import glexindex
-from .aPoly_Construction import aPoly_Construction
-from .Exploration import Exploration
from .RegressionFastARD import RegressionFastARD
from .RegressionFastLaplace import RegressionFastLaplace
from .bayes_linear import VBLinearRegression, EBLinearRegression
+warnings.filterwarnings("ignore")
+plt.style.use('../../bayesvalidrox.mplstyle')
+
class Metamodel:
-
+
def __init__(self, InputClass):
self.Inputs = InputClass
+ self.metaModel = 'PCE'
self.DisplayFlag = False
- self.slicingforValidation = False
- self.NofPa = len(InputClass.Marginals)
self.MaxPceDegree = None
self.MinPceDegree = 1
self.q = 1.0
self.P = None
self.RegMethod = None
- self.metaModel = 'PCE'
self.DimRedMethod = 'no'
- self.varPCAThreshold = 100
self.Discrepancy = None
-
self.ExpDesign = None
self.ExpDesignFlag = 'normal'
self.OptDesignFlag = False
- self.ModelOutputDict= None
+ self.ModelOutputDict = None
self.qNormDict = {}
- self.OptimalCollocationPointsBase = []
- # self.errorModel = []
- self.SeqModifiedLOO = {}
- self.seqValidError = {}
- self.SeqBME = {}
- self.SeqKLD = {}
- self.SeqDistHellinger = {}
- self.seqRMSEMean = {}
- self.seqRMSEStd = {}
- self.seqMinDist = []
self.Observations = []
self.validModelRuns = []
- # self.validSamples = []
+ self.validSamples = []
self.validLikelihoods = []
-
self.index = None
self.adaptVerbose = False
- self.Likelihoods = None
+
+ # -------------------------------------------------------------------------
+ def create_metamodel(self, Model):
+
+ self.ModelObj = Model
+ self.ExpDesign.initNrSamples = self.ExpDesign.NrSamples
+ self.NofPa = len(self.Inputs.Marginals)
+
+ if self.ExpDesign.Method == 'sequential':
+ from .sequential_design import SeqDesign
+ seq_design = SeqDesign(self)
+ metamodel = seq_design.create_PCE_SeqDesign(Model)
+
+ elif self.ExpDesign.Method == 'normal':
+ metamodel = self.create_PCE_normDesign(Model)
+
+ else:
+ raise Exception("The method for experimental design you requested"
+ " has not been implemented yet.")
+
+ # Cleanup
+ # Zip the subdirectories
+ try:
+ dir_name = Model.Name
+ key = Model.Name + '_'
+ Model.zip_subdirs(dir_name, key)
+ except:
+ pass
+
+ return metamodel
# -------------------------------------------------------------------------
class AutoVivification(dict):
@@ -119,8 +96,8 @@ class Metamodel:
except KeyError:
value = self[item] = type(self)()
return value
- # -------------------------------------------------------------------------
+ # -------------------------------------------------------------------------
def PolyBasisIndices(self, degree, q):
self.PolynomialDegrees = glexindex(start=0, stop=degree+1,
@@ -170,7 +147,7 @@ class Metamodel:
f.close()
else:
# Check if an old hdf5 file exists: if yes, rename it
- hdf5file = 'ExpDesign'+'_'+self.ModelObj.Name+'.hdf5'
+ hdf5file = 'ExpDesign'+'_'+self.ModelObj.name+'.hdf5'
if os.path.exists(hdf5file):
os.rename(hdf5file, 'old_'+hdf5file)
@@ -219,41 +196,6 @@ class Metamodel:
apolycoeffs = None
return eval_univ_basis(ED_X, n_max, polyTypes, apolycoeffs)
- # polytypes = self.ExpDesign.Polytypes
- # origSpaceDist = self.ExpDesign.JDist
- # if ED_X.ndim != 2:
- # ED_X = ED_X.reshape(1, len(ED_X))
- # NofSamples, NofPa = ED_X.shape
- # n_max = self.MaxPceDegree if n_max is None else n_max
- # P = n_max + 1
- # Allocate the output matrix
- # univ_vals = np.zeros((NofSamples, NofPa, P), dtype=np.float32)
-
- # ----------------
- # if 'arbitrary' in polytypes:
-
- # def eval_rec_rule_arbitrary(x, max_deg, polycoeffs):
- # P = max_deg+1
- # values = np.zeros((len(x), P))
-
- # for deg in range(P):
- # values[:, deg] = polyval(x, polycoeffs[deg]).T
-
- # return values
-
- # # Calculate the univariate polynomials
- # for parIdx in range(NofPa):
- # polycoeffs = self.ExpDesign.polycoeffs['p_'+str(parIdx+1)]
- # univ_vals[:, parIdx] = eval_rec_rule_arbitrary(ED_X[:, parIdx],
- # n_max,
- # polycoeffs)
- # return univ_vals
- # # ----------------
- # else:
- # from .eval_rec_rule import eval_univ_basis
- # polyTypes = self.ExpDesign.Polytypes
- # return eval_univ_basis(ED_X, n_max, polyTypes)
-
# -------------------------------------------------------------------------
def PCE_create_Psi(self, BasisIndices, univ_p_val):
"""
@@ -296,11 +238,31 @@ class Metamodel:
return Psi
# -------------------------------------------------------------------------
- def RS_Builder(self, PSI, ModelOutput, BasisIndices,
- startBasisIndices=None, RegMethod=None):
+ def RS_Builder(self, X, y, BasisIndices, RegMethod=None):
"""
+ Fit regression using the regression method provided.
+
+ Parameters
+ ----------
+ X : array-like of shape (n_samples, n_features)
+ Training vector, where n_samples is the number of samples and
+ n_features is the number of features.
+ y : array-like of shape (n_samples,)
+ Target values.
+ BasisIndices : TYPE
+ DESCRIPTION.
+ RegMethod : string, optional
+ DESCRIPTION. The default is None.
+
+ Raises
+ ------
+ ValueError
+ DESCRIPTION.
- Fit regression model.
+ Returns
+ -------
+ returnOuts : Dict
+ Fitted estimator, spareMulti-Index, sparseX and coefficients.
"""
if RegMethod is None:
@@ -313,14 +275,14 @@ class Metamodel:
compute_score = True if self.DisplayFlag else False
# inverse of the observed variance of the data
- if np.var(ModelOutput) != 0:
- Lambda = 1 / np.var(ModelOutput)
+ if np.var(y) != 0:
+ Lambda = 1 / np.var(y)
else:
Lambda = 1e-6
# Bayes sparse adaptive aPCE
if RegMethod.lower() != 'ols':
- if RegMethod.lower() == 'brr' or np.all(ModelOutput == 0):
+ if RegMethod.lower() == 'brr' or np.all(y == 0):
clf_poly = lm.BayesianRidge(n_iter=1000, tol=1e-7,
fit_intercept=True,
normalize=True,
@@ -338,8 +300,7 @@ class Metamodel:
lambda_1=Lambda, lambda_2=Lambda)
elif RegMethod.lower() == 'fastard':
- clf_poly = RegressionFastARD(start=startBasisIndices,
- fit_intercept=True,
+ clf_poly = RegressionFastARD(fit_intercept=True,
normalize=True,
compute_score=compute_score,
n_iter=300, tol=1e-10)
@@ -366,7 +327,7 @@ class Metamodel:
clf_poly = EBLinearRegression(optimizer='em')
# Fit
- clf_poly.fit(PSI, ModelOutput)
+ clf_poly.fit(X, y)
# Select the nonzero entries of coefficients
# The first column must be kept (For mean calculations)
@@ -376,269 +337,288 @@ class Metamodel:
nnz_idx = np.insert(np.nonzero(clf_poly.coef_)[0], 0, 0)
# Remove the zero entries for Bases and PSI if need be
if sparsity:
- PolynomialDegrees_Sparse = BasisIndices.toarray()[nnz_idx]
+ sparse_basis_indices = BasisIndices.toarray()[nnz_idx]
else:
- PolynomialDegrees_Sparse = BasisIndices[nnz_idx]
- PSI_Sparse = PSI[:, nnz_idx]
+ sparse_basis_indices = BasisIndices[nnz_idx]
+ sparse_X = X[:, nnz_idx]
# Store the coefficients of the regression model
- clf_poly.fit(PSI_Sparse, ModelOutput)
+ clf_poly.fit(sparse_X, y)
coeffs = clf_poly.coef_
else:
# This is for the case where all outputs are zero, thereby
# all coefficients are zero
if sparsity:
- PolynomialDegrees_Sparse = BasisIndices.toarray()
+ sparse_basis_indices = BasisIndices.toarray()
else:
- PolynomialDegrees_Sparse = BasisIndices
- PSI_Sparse = PSI
+ sparse_basis_indices = BasisIndices
+ sparse_X = X
coeffs = clf_poly.coef_
- # Raise an error, if all coeffs are zero
- # raise ValueError('All the coefficients of the regression
- # model are zero!')
# Ordinary least square method (OLS)
else:
if sparsity:
- PolynomialDegrees_Sparse = BasisIndices.toarray()
+ sparse_basis_indices = BasisIndices.toarray()
else:
- PolynomialDegrees_Sparse = BasisIndices
- PSI_Sparse = PSI
+ sparse_basis_indices = BasisIndices
+ sparse_X = X
- PsiTPsi = np.dot(PSI_Sparse.T, PSI_Sparse)
+ X_T_X = np.dot(sparse_X.T, sparse_X)
- if np.linalg.cond(PsiTPsi) > 1e-12 and \
- np.linalg.cond(PsiTPsi) < 1 / sys.float_info.epsilon:
+ if np.linalg.cond(X_T_X) > 1e-12 and \
+ np.linalg.cond(X_T_X) < 1 / sys.float_info.epsilon:
# faster
- coeffs = sp.linalg.solve(PsiTPsi, np.dot(PSI_Sparse.T, ModelOutput))
+ coeffs = sp.linalg.solve(X_T_X, np.dot(sparse_X.T, y))
else:
# stabler
- coeffs = np.dot(np.dot(np.linalg.pinv(PsiTPsi), PSI_Sparse.T),
- ModelOutput)
+ coeffs = np.dot(np.dot(np.linalg.pinv(X_T_X), sparse_X.T), y)
- return PolynomialDegrees_Sparse, PSI_Sparse, coeffs, clf_poly
+ # Create a dict to pass the outputs
+ returnOuts = dict()
+ returnOuts['clf_poly'] = clf_poly
+ returnOuts['spareMulti-Index'] = sparse_basis_indices
+ returnOuts['sparePsi'] = sparse_X
+ returnOuts['coeffs'] = coeffs
+ return returnOuts
# --------------------------------------------------------------------------------------------------------
- def adaptiveSparseaPCE(self, CollocationPoints, ModelOutput, varIdx, prevBasisIndices=None, verbose=False):
-
- NrSamples, NofPa = CollocationPoints.shape
+ def adaptiveSparseaPCE(self, ED_X, ED_Y, varIdx, verbose=False):
+ """
+ Adaptively fits the PCE model by comparing the scores of different
+ degrees and q-norm.
+
+ Parameters
+ ----------
+ ED_X : array-like of shape (n_samples, n_params)
+ Experimental design.
+ ED_Y : array-like of shape (n_samples,)
+ Target values, i.e. simulation results for the Experimental design.
+ varIdx : int
+ Index of the output.
+ verbose : bool, optional
+ Print out summary. The default is False.
+
+ Returns
+ -------
+ returnVars : Dict
+ Fitted estimator, best degree, best q-norm, LOOCVScore and
+ coefficients.
+
+ """
+
+ NrSamples, NofPa = ED_X.shape
# Initialization
qAllCoeffs, AllCoeffs = {}, {}
- qAllIndices_Sparse ,AllIndices_Sparse = {}, {}
+ qAllIndices_Sparse, AllIndices_Sparse = {}, {}
qAllclf_poly, Allclf_poly = {}, {}
qAllnTerms, AllnTerms = {}, {}
qAllLCerror, AllLCerror = {}, {}
-
+
# Extract degree array and qnorm array
- DegreeArray = np.array([*self.allBasisIndices],dtype=int)
+ DegreeArray = np.array([*self.allBasisIndices], dtype=int)
qnorm = [*self.allBasisIndices[str(int(DegreeArray[0]))]]
+ # Some options for EarlyStop
errorIncreases = False
- # Overall score indicator: this will be used to compare between different
- # methods
- bestScore = -np.inf
-
- # some options for EarlyStop
- # stop degree, if LOO error does not decrease n_checks_degree times
+ # Stop degree, if LOO error does not decrease n_checks_degree times
n_checks_degree = 3
- # stop qNorm, if criterion isn't fulfilled n_checks_qNorm times
+ # Stop qNorm, if criterion isn't fulfilled n_checks_qNorm times
n_checks_qNorm = 2
nqnorms = len(qnorm)
qNormEarlyStop = True
if nqnorms < n_checks_qNorm+1:
qNormEarlyStop = False
-
- #==============================================================================
- #==== basis adaptive polynomial chaos: repeat the calculation by
- #==== increasing polynomial degree until the highest accuracy is reached
- #==============================================================================
+
+ # =====================================================================
+ # basis adaptive polynomial chaos: repeat the calculation by increasing
+ # polynomial degree until the highest accuracy is reached
+ # =====================================================================
# For each degree check all q-norms and choose the best one
- scores = -np.inf*np.ones(DegreeArray.shape[0])
- qNormScores = -np.inf*np.ones(nqnorms)
- scores_BRR = -np.inf*np.ones(DegreeArray.shape[0])
-
+ scores = -np.inf * np.ones(DegreeArray.shape[0])
+ qNormScores = -np.inf * np.ones(nqnorms)
+
for degIdx, deg in enumerate(DegreeArray):
for qidx, q in enumerate(qnorm):
- # Extract the polynomial basis indices from the pool of allBasisIndices
+ # Extract the polynomial basis indices from the pool of
+ # allBasisIndices
BasisIndices = self.allBasisIndices[str(deg)][str(q)]
# Assemble the Psi matrix
Psi = self.PCE_create_Psi(BasisIndices, self.univ_p_val)
- # ---- Calulate the cofficients of aPCE ----------------------
- sparseBasisIndices, sparsePSI, Coeffs, clf_poly = self.RS_Builder(Psi,ModelOutput,BasisIndices)
+ # Calulate the cofficients of the meta model
+ outs = self.RS_Builder(Psi, ED_Y, BasisIndices)
# Calculate and save the score of LOOCV
- score, LCerror = self.CorrectedLeaveoutError(sparsePSI, Coeffs, ModelOutput, clf_poly)
+ score, LCerror = self.CorrectedLeaveoutError(outs['clf_poly'],
+ outs['sparePsi'],
+ outs['coeffs'],
+ ED_Y)
# Check the convergence of noise for FastARD
- if self.RegMethod == 'FastARD' and clf_poly.alpha_ < np.finfo(np.float32).eps:
+ if self.RegMethod == 'FastARD' and \
+ outs['clf_poly'].alpha_ < np.finfo(np.float32).eps:
score = -np.inf
qNormScores[qidx] = score
- qAllCoeffs[str(qidx+1)] = Coeffs
- qAllIndices_Sparse[str(qidx+1)] = sparseBasisIndices
- qAllclf_poly[str(qidx+1)] = clf_poly
+ qAllCoeffs[str(qidx+1)] = outs['coeffs']
+ qAllIndices_Sparse[str(qidx+1)] = outs['spareMulti-Index']
+ qAllclf_poly[str(qidx+1)] = outs['clf_poly']
qAllnTerms[str(qidx+1)] = BasisIndices.shape[0]
qAllLCerror[str(qidx+1)] = LCerror
# EarlyStop check
# if there are at least n_checks_qNorm entries after the
# best one, we stop
- if qNormEarlyStop and sum(np.isfinite(qNormScores)) > n_checks_qNorm:
- # If the error has increased the last two iterations, stop
- deltas = np.sign(np.diff(qNormScores[np.isfinite(qNormScores)]))
+ if qNormEarlyStop and \
+ sum(np.isfinite(qNormScores)) > n_checks_qNorm:
+ # If the error has increased the last two iterations, stop!
+ qNormScores_nonInf = qNormScores[np.isfinite(qNormScores)]
+ deltas = np.sign(np.diff(qNormScores_nonInf))
if sum(deltas[-n_checks_qNorm+1:]) == 2:
# stop the q-norm loop here
break
- if np.var(ModelOutput) == 0:
+ if np.var(ED_Y) == 0:
break
-
- # Leave the loop, if FastARD did not converge.
- # if self.RegMethod == 'FastARD' and not clf_poly.converged and deg != 1:
- # print("Degree {0} not converged!".format(deg))
- # # scores[degIdx] = -1*np.inf
- # # break
- # else:
+
# Store the score in the scores list
- bestqIdx = np.nanargmax(qNormScores)
- scores[degIdx] = qNormScores[bestqIdx] #np.max(qNormScores)
-
- AllCoeffs[str(degIdx+1)] = qAllCoeffs[str(bestqIdx+1)]
- AllIndices_Sparse[str(degIdx+1)] = qAllIndices_Sparse[str(bestqIdx+1)]
- Allclf_poly[str(degIdx+1)] = qAllclf_poly[str(bestqIdx+1)]
- AllnTerms[str(degIdx+1)] = qAllnTerms[str(bestqIdx+1)]
- AllLCerror[str(degIdx+1)] = qAllLCerror[str(bestqIdx+1)]
-
-
- # check the direction of the error (on average):
+ best_q = np.nanargmax(qNormScores)
+ scores[degIdx] = qNormScores[best_q]
+
+ AllCoeffs[str(degIdx+1)] = qAllCoeffs[str(best_q+1)]
+ AllIndices_Sparse[str(degIdx+1)] = qAllIndices_Sparse[str(best_q+1)]
+ Allclf_poly[str(degIdx+1)] = qAllclf_poly[str(best_q+1)]
+ AllnTerms[str(degIdx+1)] = qAllnTerms[str(best_q+1)]
+ AllLCerror[str(degIdx+1)] = qAllLCerror[str(best_q+1)]
+
+ # Check the direction of the error (on average):
# if it increases consistently stop the iterations
- if len(scores[scores!=-np.inf]) > n_checks_degree:
- ss = np.sign(scores[scores!=-np.inf] - np.max(scores[scores!=-np.inf])) #ss<0 error decreasing
+ if len(scores[scores != -np.inf]) > n_checks_degree:
+ scores_nonInf = scores[scores != -np.inf]
+ ss = np.sign(scores_nonInf - np.max(scores_nonInf))
+ # ss<0 error decreasing
errorIncreases = np.sum(np.sum(ss[-2:])) <= -1*n_checks_degree
-
+
if errorIncreases:
break
-
- # Check only one degree
- if np.var(ModelOutput) == 0:
+
+ # Check only one degree, if target matrix has zero variance
+ if np.var(ED_Y) == 0:
break
# ------------------ Summary of results ------------------
# Select the one with the best score and save the necessary outputs
- bestIdx = np.nanargmax(scores)+1
- coeffs = AllCoeffs[str(bestIdx)]
- PolynomialDegrees = AllIndices_Sparse[str(bestIdx)]
- clf_poly = Allclf_poly[str(bestIdx)]
+ best_deg = np.nanargmax(scores)+1
+ coeffs = AllCoeffs[str(best_deg)]
+ PolynomialDegrees = AllIndices_Sparse[str(best_deg)]
+ clf_poly = Allclf_poly[str(best_deg)]
LOOCVScore = np.nanmax(scores)
- P = AllnTerms[str(bestIdx)]
- LCerror = AllLCerror[str(bestIdx)]
+ P = AllnTerms[str(best_deg)]
+ LCerror = AllLCerror[str(best_deg)]
degree = DegreeArray[np.nanargmax(scores)]
- qnorm = float(qnorm[bestqIdx])
+ qnorm = float(qnorm[best_q])
- # ------------------ Summary of results ------------------
+ # ------------------ Print out Summary of results ------------------
if self.DisplayFlag:
# Create PSI_Sparse by removing redundent terms
nnz_idx = np.nonzero(coeffs)[0]
BasisIndices_Sparse = PolynomialDegrees[nnz_idx]
- print('Output variable %s:'%(varIdx+1))
- print('The estimation of PCE coefficients converged at polynomial degree '
- '%s with %s terms (Sparsity index = %s).'%(DegreeArray[bestIdx-1],
- len(BasisIndices_Sparse),
- round(len(BasisIndices_Sparse)/P, 3)))
- print('Final ModLOO error estimate: %.3e'%(1-max(scores)))
+ print(f'Output variable {varIdx+1}:')
+ print('The estimation of PCE coefficients converged at polynomial '
+ f'degree {DegreeArray[best_deg-1]} with '
+ f'{len(BasisIndices_Sparse)} terms (Sparsity index = '
+ f'{round(len(BasisIndices_Sparse)/P, 3)}).')
+
+ print(f'Final ModLOO error estimate: {1-max(scores):.3e}')
print('\n'+'-'*50)
-
+
if verbose:
print('='*50)
- print(' '* 10 +' Summary of results ')
+ print(' '*10 + ' Summary of results ')
print('='*50)
-
+
print("scores:\n", scores)
- print("Best score's degree:", self.DegreeArray[bestIdx-1])
+ print("Best score's degree:", self.DegreeArray[best_deg-1])
print("NO. of terms:", len(PolynomialDegrees))
print("Sparsity index:", round(len(PolynomialDegrees)/P, 3))
print("Best Indices:\n", PolynomialDegrees)
-
+
if self.RegMethod in ['BRR', 'ARD']:
fig, ax = plt.subplots(figsize=(12, 10))
plt.title("Marginal log-likelihood")
plt.plot(clf_poly.scores_, color='navy', linewidth=2)
plt.ylabel("Score")
plt.xlabel("Iterations")
- try:
- text = "$\\alpha={:.1f}$\n$\\lambda={:.3f}$\n$L={:.1f}$".format(
- clf_poly.alpha_, clf_poly.lambda_, clf_poly.scores_[-1])
- except:
- text = "$\\alpha={:.1f}$\n$\\L={:.1f}$".format(
- clf_poly.alpha_, clf_poly.scores_[-1])
-
- plt.text(0.75, 0.5, text, fontsize=18, transform = ax.transAxes)
+ if self.RegMethod.lower() == 'bbr':
+ text = f"$\\alpha={clf_poly.alpha_:.1f}$\n"
+ f"$\\lambda={clf_poly.lambda_:.3f}$\n"
+ f"$L={clf_poly.scores_[-1]:.1f}$"
+ else:
+ text = f"$\\alpha={clf_poly.alpha_:.1f}$\n$"
+ f"\\L={clf_poly.scores_[-1]:.1f}$"
+
+ plt.text(0.75, 0.5, text, fontsize=18, transform=ax.transAxes)
plt.show()
-
- print ('='*80)
-
- # Create a dict to pass the variables
+ print('='*80)
+
+ # Create a dict to pass the outputs
returnVars = dict()
-
+ returnVars['clf_poly'] = clf_poly
returnVars['degree'] = degree
returnVars['qnorm'] = qnorm
returnVars['coeffs'] = coeffs
- returnVars['PolynomialDegrees'] = PolynomialDegrees
+ returnVars['multi_indices'] = PolynomialDegrees
returnVars['LOOCVScore'] = LOOCVScore
- returnVars['clf_poly'] = clf_poly
returnVars['LCerror'] = LCerror
-
- return returnVars
-
- #--------------------------------------------------------------------------------------------------------
- def CoeffsUpdateaPCE(self, CollocationPoints, ModelOutput, BasisIndices):
-
- # Evaluate the multivariate polynomials on CollocationPoints
- Psi = self.PCE_create_Psi(BasisIndices,self.univ_p_val)
-
- # Calulate the cofficients of surrogate model
- PolynomialDegrees, PSI_Sparse, coeffs, clf_poly = self.RS_Builder(Psi, ModelOutput, BasisIndices, RegMethod='OLS') #,startIdxes=startIdxes)
-
- # Calculate and save the score of LOOCV
- LOOCVScore, LCerror = self.CorrectedLeaveoutError(PSI_Sparse, coeffs,
- ModelOutput, clf_poly)
-
- # Create a dict to pass the variables
- returnVars = dict()
-
-
+
return returnVars
-
+
# -------------------------------------------------------------------------
- def CorrectedLeaveoutError(self, TotalPsi, Coeffs, ModelOutputs, clf):
+ def CorrectedLeaveoutError(self, clf, psi, coeffs, y):
"""
- calculate the corrected LOO error for the OLS regression on regressor
- matrix PSI that generated the coefficients.
- (based on Blatman, 2009 (PhD Thesis), pg. 115-116).
+ Calculates the corrected LOO error for the OLS regression on regressor
+ matrix PSI that generated the coefficients.
+ (based on Blatman, 2009 (PhD Thesis), pg. 115-116).
This is based on the following paper:
- ""Blatman, G., & Sudret, B. (2011). Adaptive sparse polynomial
- chaos expansion based on least angle regression.
- Journal of Computational Physics, 230(6), 2345-2367.""
+ ""Blatman, G., & Sudret, B. (2011). Adaptive sparse polynomial
+ chaos expansion based on least angle regression.
+ Journal of Computational Physics, 230(6), 2345-2367.""
+
+ Parameters
+ ----------
+ clf : object
+ Fitted estimator.
+ psi : array-like of shape (n_samples, n_features)
+ The multivariate orthogonal polynomials (regressor).
+ coeffs : array-like of shape (n_features,)
+ Estimated cofficients.
+ y : array-like of shape (n_samples,)
+ Target values.
+
+ Returns
+ -------
+ Q_2 : float
+ LOOCV Validation score (1-LOOCV erro).
+ residual : array-like of shape (n_samples,)
+ Residual values (y - predicted targets).
- Psi: the orthogonal polynomials of the response surface
"""
- Psi = np.array(TotalPsi, dtype=float)
+ psi = np.array(psi, dtype=float)
# Create PSI_Sparse by removing redundent terms
- nnz_idx = np.nonzero(Coeffs)[0]
+ nnz_idx = np.nonzero(coeffs)[0]
if len(nnz_idx) == 0:
nnz_idx = [0]
- PSI_Sparse = Psi[:, nnz_idx]
+ PSI_Sparse = psi[:, nnz_idx]
# NrCoeffs of aPCEs
P = len(nnz_idx)
# NrEvaluation (Size of experimental design)
- N = Psi.shape[0]
+ N = psi.shape[0]
# Build the projection matrix
PsiTPsi = np.dot(PSI_Sparse.T, PSI_Sparse)
@@ -661,17 +641,17 @@ class Metamodel:
# ------ Calculate Error Loocv for each measurement point ----
# Residuals
if isinstance(clf, list):
- residual = np.dot(Psi, Coeffs) - ModelOutputs
+ residual = np.dot(psi, coeffs) - y
else:
- residual = clf.predict(Psi) - ModelOutputs
+ residual = clf.predict(psi) - y
# Variance
- varY = np.var(ModelOutputs)
+ varY = np.var(y)
if varY == 0:
normEmpErr = 0
ErrLoo = 0
- LCerror = np.zeros((ModelOutputs.shape))
+ LCerror = np.zeros((y.shape))
else:
normEmpErr = np.mean(residual**2)/varY
@@ -686,7 +666,7 @@ class Metamodel:
# Corrected Error for over-determined system
trM = np.trace(M)
if trM < 0 or abs(trM) > 1e6:
- trM = np.trace(np.linalg.pinv(np.dot(Psi.T, Psi)))
+ trM = np.trace(np.linalg.pinv(np.dot(psi.T, psi)))
# Over-determined system of Equation
if N > P:
@@ -702,67 +682,118 @@ class Metamodel:
return Q_2, residual
- #--------------------------------------------------------------------------------------------------------
+ # -------------------------------------------------------------------------
def PCATransformation(self, Output):
-
- # Transform via Principal Component Analysis
- from sklearn.decomposition import PCA as sklearnPCA
- # from sklearn.decomposition import TruncatedSVD as sklearnPCA
+ """
+ Transforms the targets (outputs) via Principal Component Analysis
+
+ Parameters
+ ----------
+ Output : array-like of shape (n_samples,)
+ Target values.
+
+ Returns
+ -------
+ pca : object
+ Fitted sklearnPCA object.
+ OutputMatrix : array-like of shape (n_samples,)
+ Transformed target values.
+ """
+ # Transform via Principal Component Analysis
+ if isinstance(self, 'varPCAThreshold'):
+ varPCAThreshold = self.varPCAThreshold
+ else:
+ varPCAThreshold = 100.0
n_samples, n_features = Output.shape
- covar_matrix = sklearnPCA(n_components=None) #, svd_solver='full')
+
+ # Instantiate and fit sklearnPCA object
+ covar_matrix = sklearnPCA(n_components=None)
covar_matrix.fit(Output)
- var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_, decimals=5)*100)
+ var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_,
+ decimals=5)*100)
+ # Find the number of components to explain self.varPCAThreshold of
+ # variance
try:
- selected_n_components = np.where(var>=self.varPCAThreshold)[0][0] + 1
- except:
+ selected_n_components = np.where(var >= varPCAThreshold)[0][0] + 1
+ except ValueError:
selected_n_components = min(n_samples, n_features)
-
+
nComponents = min(n_samples, n_features, selected_n_components)
- pca = sklearnPCA(n_components=nComponents, svd_solver='randomized') #, svd_solver='full')
+
+ # Fit and transform with the selected number of components
+ pca = sklearnPCA(n_components=nComponents, svd_solver='randomized')
OutputMatrix = pca.fit_transform(Output)
-
+
return pca, OutputMatrix
- #--------------------------------------------------------------------------------------------------------
- def GaussianProcessEmulator(self, X, y, nugTerm=None,autoSelect=False, varIdx=None):
-
+ # -------------------------------------------------------------------------
+ def GaussianProcessEmulator(self, X, y, nugTerm=None, autoSelect=False,
+ varIdx=None):
+ """
+ Fits a Gaussian Process Emulator to the target given the training
+ points.
+
+ Parameters
+ ----------
+ X : array-like of shape (n_samples, n_params)
+ Training points.
+ y : array-like of shape (n_samples,)
+ Target values.
+ nugTerm : float, optional
+ Nugget term. The default is None, i.e. variance of y.
+ autoSelect : bool, optional
+ Loop over some kernels and select the best. The default is False.
+ varIdx : int, optional
+ The index number. The default is None.
+
+ Returns
+ -------
+ gp : object
+ Fitted estimator.
+
+ """
+
nugTerm = nugTerm if nugTerm else np.var(y)
- 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-15, 1e5),
+ Kernels = [nugTerm * kernels.RBF(length_scale=1.0,
+ length_scale_bounds=(1e-25, 1e15)),
+ nugTerm * kernels.RationalQuadratic(length_scale=0.2,
+ alpha=1.0),
+ nugTerm * kernels.Matern(length_scale=1.0,
+ length_scale_bounds=(1e-15, 1e5),
nu=1.5)]
- if autoSelect:# Automatic selection of the kernel
+ # Automatic selection of the kernel
+ if autoSelect:
gp = {}
BME = []
for i, kernel in enumerate(Kernels):
- gp[i] = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=2,
- normalize_y=True)
-
- # Fit to data using Maximum Likelihood Estimation of the parameters
+ gp[i] = GaussianProcessRegressor(kernel=kernel,
+ n_restarts_optimizer=3,
+ normalize_y=False)
+
+ # Fit to data using Maximum Likelihood Estimation
gp[i].fit(X, y)
-
+
+ # Store the MLE as BME score
BME.append(gp[i].log_marginal_likelihood())
-
- gp = gp[np.argmax(BME)]
-
+
+ gp = gp[np.argmax(BME)]
+
else:
- gp = GaussianProcessRegressor(kernel=Kernels[0], n_restarts_optimizer=3,
- normalize_y=False)
+ gp = GaussianProcessRegressor(kernel=Kernels[0],
+ n_restarts_optimizer=3,
+ normalize_y=False)
gp.fit(X, y)
# Compute score
if varIdx is not None:
Score = gp.score(X, y)
-
print('-'*50)
- print('Output variable %s:'%(varIdx))
+ print(f'Output variable {varIdx}:')
print('The estimation of GPE coefficients converged,')
- print('with the R^2 score: {0:.3f}'.format(Score))
- #print("BME:", gp.log_marginal_likelihood())
+ print(f'with the R^2 score: {Score:.3f}')
print('-'*50)
return gp
@@ -773,6 +804,18 @@ class Metamodel:
This function loops over the outputs and each time step/point and fits
the meta model.
+ Parameters
+ ----------
+ Model : object
+ Model object.
+ OptDesignFlag : bool, optional
+ Flag for a sequential design in silent mode. The default is False.
+
+ Returns
+ -------
+ self: object
+ Meta-model object.
+
"""
# Get the collocation points to run the forward model
@@ -791,33 +834,24 @@ class Metamodel:
self.LCerror = self.AutoVivification()
self.x_scaler = {}
+ # Read observations or MCReference
if self.ExpDesignFlag != 'sequential':
- # Read observations or MCReference
- if Model.MeasurementFile is not None \
- or len(Model.Observations) != 0:
+ if len(Model.Observations) != 0:
self.Observations = Model.read_Observation()
# Define the DegreeArray
- maxDeg = self.MaxPceDegree
nSamples, ndim = CollocationPoints.shape
- nitr = nSamples - self.ExpDesign.initNrSamples
- M_uptoMax = lambda maxDeg: np.array([math.factorial(ndim+d)//(math.factorial(ndim)*math.factorial(d)) for d in range(1,maxDeg+1)])
- d = nitr if nitr != 0 and self.NofPa > 5 else 1
- degNew = range(1,maxDeg+1)[np.argmin(abs(M_uptoMax(maxDeg)-ndim*nSamples*d))]
- degNew = degNew if self.ExpDesignFlag == 'sequential' else self.MaxPceDegree
- self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
- if degNew > self.MinPceDegree and self.RegMethod.lower() != 'fastard':
- self.DegreeArray = np.arange(self.MinPceDegree, degNew+1)
- else:
- self.DegreeArray = np.array([degNew])
+ self.DegreeArray = self.select_degree(ndim, nSamples)
# Generate all basis indices
self.allBasisIndices = self.AutoVivification()
for deg in self.DegreeArray:
- if deg not in np.fromiter(self.allBasisIndices.keys(), dtype=float):
+ keys = self.allBasisIndices.keys()
+ if deg not in np.fromiter(keys, dtype=float):
# Generate the polynomial basis indices
for qidx, q in enumerate(self.q):
- self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree=deg, q=q)
+ basis_indices = self.PolyBasisIndices(degree=deg, q=q)
+ self.allBasisIndices[str(deg)][str(q)] = basis_indices
# Evaluate the univariate polynomials on ExpDesign
if self.metaModel.lower() != 'gpe':
@@ -828,7 +862,8 @@ class Metamodel:
# --- Loop through data points and fit the surrogate ---
if not OptDesignFlag:
- print(f"\n>>>> Training the {self.metaModel} metamodel started. <<<<<<\n")
+ print(f"\n>>>> Training the {self.metaModel} metamodel "
+ "started. <<<<<<\n")
items = tqdm(OutputDict.items(), desc="Fitting regression")
else:
items = OutputDict.items()
@@ -838,9 +873,9 @@ class Metamodel:
# Dimensionality reduction with PCA, if specified
if self.DimRedMethod.lower() == 'pca':
- self.pca[key], OutputMatrix = self.PCATransformation(Output)
+ self.pca[key], target = self.PCATransformation(Output)
else:
- OutputMatrix = Output
+ target = Output
# Parallel fit regression
if self.metaModel.lower() == 'gpe':
@@ -850,38 +885,28 @@ class Metamodel:
self.x_scaler[key] = scaler
- outDict = Parallel(n_jobs=-1, prefer='threads')(
- delayed(self.GaussianProcessEmulator)(X_S, OutputMatrix[:, idx])
- for idx in range(OutputMatrix.shape[1]))
+ out = Parallel(n_jobs=-1, prefer='threads')(
+ delayed(self.GaussianProcessEmulator)(X_S, target[:, idx])
+ for idx in range(target.shape[1]))
- for idx in range(OutputMatrix.shape[1]):
- self.gp_poly[key][f"y_{idx+1}"] = outDict[idx]
+ for idx in range(target.shape[1]):
+ self.gp_poly[key][f"y_{idx+1}"] = out[idx]
else:
- outDict = Parallel(n_jobs=-1, prefer='threads')(
+ out = Parallel(n_jobs=-1, prefer='threads')(
delayed(self.adaptiveSparseaPCE)(CollocationPoints,
- OutputMatrix[:, idx], idx)
- for idx in range(OutputMatrix.shape[1]))
+ target[:, idx], idx)
+ for idx in range(target.shape[1]))
- for idx in range(OutputMatrix.shape[1]):
+ for i in range(target.shape[1]):
# Create a dict to pass the variables
- self.DegreeDict[key][f"y_{idx+1}"] = outDict[idx]['degree']
- self.qNormDict[key][f"y_{idx+1}"] = outDict[idx]['qnorm']
- self.CoeffsDict[key][f"y_{idx+1}"] = outDict[idx]['coeffs']
- self.BasisDict[key][f"y_{idx+1}"] = outDict[idx]['PolynomialDegrees']
- self.ScoresDict[key][f"y_{idx+1}"] = outDict[idx]['LOOCVScore']
- self.clf_poly[key][f"y_{idx+1}"] = outDict[idx]['clf_poly']
- self.LCerror[key][f"y_{idx+1}"] = outDict[idx]['LCerror']
-
- # PCEKriging (PCE with Gaussian Process Emulator)
- # if self.metaModel.lower() == 'pcekriging':
- # # outDict = Parallel(n_jobs=-1)(
- # # delayed(self.GaussianProcessEmulator)(CollocationPoints,self.LCerror[key]["y_"+str(idx+1)])
- # # for idx in range(OutputMatrix.shape[1]))
-
- # for idx in range(OutputMatrix.shape[1]):
- # self.gp_poly[key]["y_"+str(idx+1)] = self.GaussianProcessEmulator(CollocationPoints, self.LCerror[key]["y_"+str(idx+1)])
- # self.gp_poly[key]["y_"+str(idx+1)] = outDict[idx]
+ self.DegreeDict[key][f"y_{i+1}"] = out[i]['degree']
+ self.qNormDict[key][f"y_{i+1}"] = out[i]['qnorm']
+ self.CoeffsDict[key][f"y_{i+1}"] = out[i]['coeffs']
+ self.BasisDict[key][f"y_{i+1}"] = out[i]['multi_indices']
+ self.ScoresDict[key][f"y_{i+1}"] = out[i]['LOOCVScore']
+ self.clf_poly[key][f"y_{i+1}"] = out[i]['clf_poly']
+ self.LCerror[key][f"y_{i+1}"] = out[i]['LCerror']
if not OptDesignFlag:
print(f"\n>>>> Training the {self.metaModel} metamodel"
@@ -889,1038 +914,130 @@ class Metamodel:
return self
- #--------------------------------------------------------------------------------------------------------
- def create_ModelError(self,X,y,Name='Calib'):
+ # -------------------------------------------------------------------------
+ def select_degree(self, ndim, nSamples):
+ # Define the DegreeArray
+ maxDeg = self.MaxPceDegree
+ nitr = nSamples - self.ExpDesign.initNrSamples
+
+ # Check q-norm
+ if not np.isscalar(self.q):
+ self.q = np.array(self.q)
+ else:
+ self.q = np.array([self.q])
+
+ def M_uptoMax(maxDeg):
+ n_combo = np.zeros(maxDeg)
+ for i, d in enumerate(range(1, maxDeg+1)):
+ n_combo[i] = math.factorial(ndim+d)
+ n_combo[i] /= math.factorial(ndim) * math.factorial(d)
+ return n_combo
+
+ if self.ExpDesignFlag != 'sequential':
+ degNew = self.MaxPceDegree
+ else:
+ d = nitr if nitr != 0 and self.NofPa > 5 else 1
+ min_index = np.argmin(abs(M_uptoMax(maxDeg)-ndim*nSamples*d))
+ degNew = range(1, maxDeg+1)[min_index]
+
+ if degNew > self.MinPceDegree and self.RegMethod.lower() != 'fastard':
+ DegreeArray = np.arange(self.MinPceDegree, degNew+1)
+ else:
+ DegreeArray = np.array([degNew])
+
+ return DegreeArray
+
+ # -------------------------------------------------------------------------
+ def create_ModelError(self, X, y, Name='Calib'):
"""
- This function loops over the outputs and each time step/point to compute the PCE
- coefficients.
-
+ This function fits a GPE-based model error.
+
+ Parameters
+ ----------
+ X : array-like of shape (n_outputs, n_inputs)
+ Input array. It can contain any forcing inputs or coordinates of
+ extracted data.
+ y : array-like of shape (n_outputs,)
+ The model response for the MAP parameter set.
+ Name : string, optional
+ Calibration or validation. The default is 'Calib'.
+
+ Returns
+ -------
+ TYPE
+ DESCRIPTION.
+
"""
Model = self.ModelObj
outputNames = Model.Output.Names
self.errorRegMethod = 'GPE'
self.errorclf_poly = self.AutoVivification()
self.errorScale = self.AutoVivification()
-
+
# Read data
if Name.lower() == 'calib':
MeasuredData = Model.read_Observation()
else:
MeasuredData = Model.read_ObservationValid()
-
+
# Fitting GPR based bias model
for out in outputNames:
+ nan_idx = ~np.isnan(MeasuredData[out])
# Select data
try:
- data = MeasuredData[out].to_numpy()[~np.isnan(MeasuredData[out])]
+ data = MeasuredData[out].values[nan_idx]
except:
- data = MeasuredData[out][~np.isnan(MeasuredData[out])]
-
+ data = MeasuredData[out][nan_idx]
+
# Prepare the input matrix
scaler = MinMaxScaler()
- delta = data #- y[out][0]
- BiasInputs = np.hstack((X[out],y[out].reshape(-1,1)))
+ delta = data # - y[out][0]
+ BiasInputs = np.hstack((X[out], y[out].reshape(-1, 1)))
X_S = scaler.fit_transform(BiasInputs)
-
+ gp = self.GaussianProcessEmulator(X_S, delta)
+
self.errorScale[out]["y_1"] = scaler
- self.errorclf_poly[out]["y_1"] = self.GaussianProcessEmulator(X_S, delta,nugTerm=np.var(delta))
-
+ self.errorclf_poly[out]["y_1"] = gp
+
return self
-
- #--------------------------------------------------------------------------------------------------------
- def eval_errormodel(self,X,y_pred):
+
+ # -------------------------------------------------------------------------
+ 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])
PCEOutputs_std = np.zeros_like(y_pred[Outkey])
-
+
for Inkey, InIdxValues in ValuesDict.items():
-
+
gp = self.errorclf_poly[Outkey][Inkey]
scaler = self.errorScale[Outkey][Inkey]
-
+
# Transform Samples using scaler
- for j,pred in enumerate(y_pred[Outkey]):
- BiasInputs = np.hstack((X[Outkey],pred.reshape(-1,1)))
+ for j, pred in enumerate(y_pred[Outkey]):
+ BiasInputs = 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
+ y_hat, y_std = gp.predict(Samples_S, return_std=True)
+ PCEOutputs_mean[j] = y_hat
+ PCEOutputs_std[j] = y_std
+ # PCEOutputs_mean[j] += pred
meanPCEOutputs[Outkey] = PCEOutputs_mean
stdPCEOutputs[Outkey] = PCEOutputs_std
-
- return meanPCEOutputs, stdPCEOutputs
-
- #--------------------------------------------------------------------------------------------------------
- def util_VarBasedDesign(self, X_can, index, UtilMethod='Entropy'):
- """
- Computes the exploitation scores based on:
- active learning MacKay(ALM) and active learning Cohn (ALC)
- Paper: Sequential Design with Mutual Information for Computer Experiments (MICE):
- Emulation of a Tsunami Model by Beck and Guillas (2016)
- """
- OutputDictY = self.ExpDesign.Y
- OutputNames = list(OutputDictY.keys())[1:]
-
- if UtilMethod == 'Entropy':
- # ----- Entropy/MMSE/active learning MacKay(ALM) -----
- # Compute perdiction variance of the old model
- Y_PC_can, std_PC_can = self.eval_metamodel(samples=X_can)
- canPredVar = {key:std_PC_can[key]**2 for key in OutputNames}
-
- varPCE = np.zeros((len(OutputNames),X_can.shape[0]))
- for KeyIdx, key in enumerate(OutputNames):
- varPCE[KeyIdx] = np.max(canPredVar[key],axis=1)
- score = np.max(varPCE, axis=0)
-
- elif UtilMethod == 'EIGF':
- # ----- Expected Improvement for Global fit -----
- # Eq (5) from Liu et al.(2018)
- # Compute perdiction error and variance of the old model
- 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}
-
- EIGF_PCE = np.zeros((len(OutputNames),X_can.shape[0]))
- for KeyIdx, key in enumerate(OutputNames):
- residual = predError[key] - OutputDictY[key][index]
- var = canPredVar[key]
- EIGF_PCE[KeyIdx] = np.max(residual**2 + var,axis=1)
- score = np.max(EIGF_PCE, axis=0)
-
- return -1 * score # -1 is for minimization instead of maximization
-
- #--------------------------------------------------------------------------------------------------------
- def util_BayesianActiveDesign(self, X_can, Model, sigma2Dict, var='DKL'):
-
- nTriningPoints = self.ExpDesign.X.shape[0]
-
- # Evaluate the PCE metamodels at that location ???
- Y_mean_can, Y_std_can = self.eval_metamodel(samples=np.array([X_can]))
-
- # Get the data
- ObservationData = self.Observations
- nObs = Model.nObs
- # TODO: Analytical DKL
- # 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_PC_MC, std_PC_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_PC_MC[key] = np.zeros((MCsize,nObs))
- logsamples = np.zeros((MCsize,nObs))
- for i in range(nObs):
- NormalDensity = stats.norm(means[i], stds[i])
- Y_PC_MC[key][:,i] = NormalDensity.rvs(MCsize)
- logsamples[:,i] = NormalDensity.logpdf(Y_PC_MC[key][:,i])
-
- logPriorLikelihoods = np.sum(logsamples, axis=1)
- std_PC_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_PC_MC,std_PC_MC, ObservationData, 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
-
- # Rejection Step
- # Random numbers between 0 and 1
- unif = np.random.rand(1, MCsize)[0]
-
- # Reject the poorly performed prior
- accepted = (Likelihoods/np.max(Likelihoods)) >= unif
-
- # Prior-based estimation of BME
- logBME = np.log(np.nanmean(Likelihoods))
-
- # Posterior-based expectation of likelihoods
- postLikelihoods = Likelihoods[accepted]
- postExpLikelihoods = np.mean(np.log(postLikelihoods))
-
- # Posterior-based expectation of prior densities
- postExpPrior = np.mean(logPriorLikelihoods[accepted])
-
- # Utility function Eq.2 in Ref. (2)
- # Posterior covariance matrix after observing data y
- # Kullback-Leibler Divergence (Sergey's paper)
- if var == 'DKL':
-
- # TODO: Calculate the correction factor for BME
-# BMECorrFactor = self.BME_Corr_Weight(PCE_SparseBayes_can, ObservationData, sigma2Dict)
-# BME += BMECorrFactor
- # Haun et al implementation
- # U_J_d = np.mean(np.log(Likelihoods[Likelihoods!=0])- logBME)
- U_J_d = postExpLikelihoods - logBME
-
-
- # Marginal log likelihood
- elif var == 'BME':
- U_J_d = logBME
-
- # Entropy-based information gain
- elif var == 'infEntropy':
- logBME = np.log(np.nanmean(Likelihoods))
- infEntropy = logBME - postExpPrior - postExpLikelihoods
- U_J_d = infEntropy * -1 # -1 for minimization
-
- # Bayesian information criterion
- elif var == 'BIC':
- nModelParams = max(len(v) for val in self.CoeffsDict.values() for v in val.values())
- maxL = np.nanmax(Likelihoods)
- U_J_d = -2* np.log(maxL) + np.log(nObs) * nModelParams
-
- # Akaike information criterion
- elif var == 'AIC':
- nModelParams = max(len(v) for val in self.CoeffsDict.values() for v in val.values())
- maxlogL = np.log(np.nanmax(Likelihoods))
- AIC = -2 * maxlogL + 2 * nModelParams
- penTerm = 0 #2 * nModelParams * (nModelParams+1) / (nObs-nModelParams-1)
- U_J_d = 1*(AIC + penTerm)
-
- # Deviance information criterion
- elif var == 'DIC':
- # 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, ObservationData,
- sigma2Dict)
- DIC = -2 * np.log(Likelihoods_theta_mean) + 2 * N_star_p
-
- U_J_d = DIC
-
- else:
- print('The algorithm you requested has not been implemented yet!')
-
- # Handle inf and NaN (replace by zero)
- if np.isnan(U_J_d) or U_J_d==-np.inf or U_J_d==np.inf: U_J_d = 0.0
-
- return -1 * U_J_d # -1 is for minimization instead of maximization
-
- #--------------------------------------------------------------------------------------------------------
- def util_BayesianDesign(self,X_can, X_MC, Model, sigma2Dict, var='DKL'):
-
- # To avoid changes ub original aPCE object
- from copy import deepcopy
- PCEModel = deepcopy(self)
-
- # Old Experimental design
- oldExpDesignX = PCEModel.ExpDesign.X
- oldExpDesignY = PCEModel.ExpDesign.Y
-
- # Evaluate the PCE metamodels at that location ???
- Y_PC_can, _ = self.eval_metamodel(samples=np.array([X_can]))
-
- # Add all suggestion as new 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]
-
- PCEModel.ExpDesign.SamplingMethod = 'user'
- PCEModel.ExpDesign.X = NewExpDesignX
- PCEModel.ExpDesign.Y = NewExpDesignY
-
- # Create the SparseBayes-based PCE metamodel:
- PCEModel.Inputs.polycoeffsFlag = False
- univ_p_val = self.univ_basis_vals(X_can)
- G_n_m_all = np.zeros((len(Model.Output.Names),Model.nObs))
-
- for idx, key in enumerate(Model.Output.Names):
- for i in range(Model.nObs):
- BasisIndices = PCEModel.BasisDict[key]["y_"+str(i+1)]
- clf_poly = PCEModel.clf_poly[key]["y_"+str(i+1)]
- Mn = clf_poly.coef_
- Sn = clf_poly.sigma_
- beta = clf_poly.alpha_
- Psi = self.PCE_create_Psi(BasisIndices,univ_p_val)
- Sn_new_inv = np.linalg.inv(Sn) + beta * np.dot(Psi[:,clf_poly.active_].T,Psi[:,clf_poly.active_])
- Sn_new = np.linalg.inv(Sn_new_inv)
-
- Mn_new = np.dot(Sn_new_inv, Mn[clf_poly.active_]).reshape(-1,1)
- Mn_new += beta * np.dot(Psi[:,clf_poly.active_].T, Y_PC_can[key][0,i])
- Mn_new = np.dot(Sn_new,Mn_new).flatten()
-
- # Compute new moments
- mean_old = Mn[0]
- mean_new = Mn_new[0]
- std_old = np.sqrt(np.sum(np.square(Mn[1:])))
- std_new = np.sqrt(np.sum(np.square(Mn_new[1:])))
-
- G_n_m = np.log(std_old/std_new) - 1./2
- G_n_m += std_new**2 / (2*std_new**2)
- G_n_m += (mean_new - mean_old)**2 / (2*std_old**2)
-
- G_n_m_all[idx,i] = G_n_m
-
-
- clf_poly.coef_[clf_poly.active_] = Mn_new
- clf_poly.sigma_ = Sn_new
- PCEModel.clf_poly[key]["y_"+str(i+1)] = clf_poly
-
- # return np.sum(G_n_m_all)
- PCE_SparseBayes_can = PCEModel#.create_PCE_normDesign(Model, OptDesignFlag=True)
-
- # Get the data
- ObservationData = self.Observations
-
- # ---------- Inner MC simulation for computing Utility Value ----------
- # Estimation of the integral via Monte Varlo integration
- MCsize = X_MC.shape[0]
- ESS = 0
-
- while ((ESS > MCsize) or (ESS < 1)):
-
- # Enriching Monte Carlo samples if need be
- if ESS != 0:
- X_MC = self.ExpDesign.GetSample(MCsize, 'random');
-
- # Evaluate the PCEModel at the given samples
- 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)
- Likelihoods = self.normpdf(Y_PC_MC,std_PC_MC, ObservationData, sigma2Dict)
-
- # Check the Effective Sample Size (1<ESS<MCsize)
- ESS = 1 / np.sum(np.square(Likelihoods/np.sum(Likelihoods)))
-
- # Enlarge sample size if it doesn't fulfill the criteria
- if ((ESS > MCsize) or (ESS < 1)):
- MCsize *= 10
- ESS = 0
-
- # Rejection Step
- # Random numbers between 0 and 1
- unif = np.random.rand(1, MCsize)[0]
-
- # Reject the poorly performed prior
- accepted = (Likelihoods/np.max(Likelihoods)) >= unif
-
-
- # -------------------- Utility functions --------------------
- # Utility function Eq.2 in Ref. (2)
- # Kullback-Leibler Divergence (Sergey's paper)
- if var == 'DKL':
-
- # Prior-based estimation of BME
- logBME = np.log(np.nanmean(Likelihoods))
-
- # Posterior-based expectation of likelihoods
- postLikelihoods = Likelihoods[accepted]
- postExpLikelihoods = np.mean(np.log(postLikelihoods))
-
- # Haun et al implementation
- U_J_d = np.mean(np.log(Likelihoods[Likelihoods!=0])- logBME)
-
- U_J_d = np.sum(G_n_m_all)
- # Ryan et al (2014) implementation
- # importanceWeights = Likelihoods[Likelihoods!=0]/np.sum(Likelihoods[Likelihoods!=0])
- # U_J_d = np.mean(importanceWeights*np.log(Likelihoods[Likelihoods!=0])) - logBME
-
- # U_J_d = postExpLikelihoods - logBME
-
- # Marginal likelihood
- elif var == 'BME':
-
- # Prior-based estimation of BME
- logBME = np.log(np.nanmean(Likelihoods))
-
- U_J_d = logBME
-
- # Bayes risk likelihood
- elif var == 'BayesRisk':
-
- U_J_d = -1 * np.var(Likelihoods)
-
- # Entropy-based information gain
- elif var == 'infEntropy':
- # Prior-based estimation of BME
- logBME = np.log(np.nanmean(Likelihoods))
-
- # Posterior-based expectation of likelihoods
- postLikelihoods = Likelihoods[accepted] / np.nansum(Likelihoods[accepted])
- postExpLikelihoods = np.mean(np.log(postLikelihoods))
-
- # Posterior-based expectation of prior densities
- postExpPrior = np.mean(logPriorLikelihoods[accepted])
-
-
- infEntropy = logBME - postExpPrior - postExpLikelihoods
-
- U_J_d = infEntropy * -1 # -1 for minimization
-
- # D-Posterior-precision
- elif var == 'DPP':
- X_Posterior = X_MC[accepted]
- # covariance of the posterior parameters
- U_J_d = -np.log(np.linalg.det(np.cov(X_Posterior)))
-
- # A-Posterior-precision
- elif var == 'APP':
- X_Posterior = X_MC[accepted]
- # trace of the posterior parameters
- U_J_d = -np.log(np.trace(np.cov(X_Posterior)))
-
- else:
- print('The algorithm you requested has not been implemented yet!')
-
- return -1 * U_J_d # -1 is for minimization instead of maximization
-
- #--------------------------------------------------------------------------------------------------------
- def subdomain(self, Bounds, NrofNewSample):
-
- NofPa = self.NofPa
- NofSubdomain = NrofNewSample+1
- LinSpace = np.zeros((NofPa,NofSubdomain))
-
- for i in range(NofPa):
- LinSpace[i]=np.linspace(start = Bounds[i][0], stop = Bounds[i][1], num = NofSubdomain)
- Subdomains=[]
-
- for k in range(NofSubdomain-1):
- mylist = []
- for i in range(NofPa):
- mylist.append((LinSpace[i,k+0],LinSpace[i,k+1]))
- Subdomains.append(tuple(mylist))
- return Subdomains
-
- #--------------------------------------------------------------------------------------------------------
- def Utility_runner(self, method, Model, allCandidates,index, sigma2Dict=None, var=None, X_MC=None):
-
- if method == 'VarOptDesign':
- U_J_d = self.util_VarBasedDesign(allCandidates, index, var)
-
- elif method == 'BayesActDesign':
- NCandidate = allCandidates.shape[0]
- U_J_d = np.zeros((NCandidate))
- for idx, X_can in tqdm(enumerate(allCandidates), ascii=True, desc ="OptBayesianDesign"):
- U_J_d[idx] = self.util_BayesianActiveDesign(X_can, Model, sigma2Dict, var)
-
- elif method == 'BayesOptDesign':
- NCandidate = allCandidates.shape[0]
- U_J_d = np.zeros((NCandidate))
- for idx, X_can in tqdm(enumerate(allCandidates), ascii=True, desc ="OptBayesianDesign"):
- U_J_d[idx] = self.util_BayesianDesign(X_can,X_MC, Model, sigma2Dict, var)
-
- return (index,-1 * U_J_d)
- #--------------------------------------------------------------------------------------------------------
-
- def dual_annealing(self,method,Bounds,sigma2Dict, var, Run_No, debugFlag=False):
-
- Model = self.ModelObj
- MaxFunItr = self.ExpDesign.MaxFunItr
-
- import scipy.optimize as opt
- if method == 'VarOptDesign':
- Res_Global = opt.dual_annealing(self.util_VarBasedDesign, bounds=Bounds,
- args=(Model,var,), maxfun=MaxFunItr)
-
- elif method == 'BayesOptDesign':
- Res_Global = opt.dual_annealing(self.util_BayesianDesign, bounds=Bounds,
- args=(Model,sigma2Dict,var), maxfun=MaxFunItr)
-
- if debugFlag:
- print("global minimum: xmin = {0}, f(xmin) = {1:.6f}, nfev = {2}".format(
- Res_Global.x, Res_Global.fun, Res_Global.nfev))
-
- return (Run_No,Res_Global.x)
- #--------------------------------------------------------------------------------------------------------
- def tradoffWeights(self, TradeOffScheme, OldExpDesign, OutputDictY):
-
- if TradeOffScheme == 'None':
- weightExploration = 0
-
- elif TradeOffScheme == 'equal':
- weightExploration = 0.5
-
- elif TradeOffScheme == 'epsilon-decreasing':
- # epsilon-decreasing scheme
- # Start with more exploration and increase the influence of exploitation
- # along the way with a exponential decay function
- initNSamples, maxNSamples = self.ExpDesign.initNrSamples, self.ExpDesign.MaxNSamples
- itrNumber = (self.ExpDesign.X.shape[0] - initNSamples)//self.ExpDesign.NrofNewSample
-
- from scipy import signal
- tau2 = -(maxNSamples-initNSamples-1) / np.log(1e-8)
- weightExploration = signal.exponential(maxNSamples-initNSamples, 0, tau2, False)[itrNumber]
-
-
- elif TradeOffScheme == 'adaptive':
-
- # Extract itrNumber
- initNSamples, maxNSamples = self.ExpDesign.initNrSamples, self.ExpDesign.MaxNSamples
- itrNumber = (self.ExpDesign.X.shape[0] - initNSamples)//self.ExpDesign.NrofNewSample
-
- if itrNumber == 0:
- weightExploration = 0.5
- else:
- # # Extract the model errors from the last and next to last iterations
- # errorModel_i , errorModel_i_1 = self.errorModel[itrNumber], self.errorModel[itrNumber-1]
-
- # # Evaluate the error models for all selected samples so far
- # eLCAllCands_i, _ = errorModel_i.eval_errormodel(OldExpDesign)
- # eLCAllCands_i_1, _ = errorModel_i_1.eval_errormodel(OldExpDesign)
-
-
- # # Local improvement of LC error at last selected design
- # sl_i = np.max(np.dstack(eLCAllCands_i.values())[-1])
- # sl_i_1 = np.max(np.dstack(eLCAllCands_i_1.values())[-1])
-
- # p = sl_i**2 / sl_i_1**2
-
- # # Global improvement of LC error at OldExpDesign
- # sg_i = np.max(np.dstack(eLCAllCands_i.values()),axis=1)
- # sg_i_1 = np.max(np.dstack(eLCAllCands_i_1.values()),axis=1)
-
- # q = np.sum(np.square(sg_i)) / np.sum(np.square(sg_i_1))
-
- # weightExploration = min([0.5*p/q, 1])
-
- # TODO: New adaptive trade-off according to Liu et al. (2017)
- # Mean squared error for last design point
- lastPCEY, _ = self.eval_metamodel(samples=OldExpDesign[-1].reshape(1,-1))
-
- mseError = mean_squared_error(np.array(list(lastPCEY.values()))[:,0], np.array(list(OutputDictY.values())[1:])[:,-1,:])
-
- # Mean squared CV - error for last design point
- mseCVError = np.mean(np.square(([v[-1] for V in self.LCerror.values() for v in V.values()])))
-
- weightExploration = 0.99 * min([0.5*mseError/mseCVError, 1])
-
- return weightExploration
- #--------------------------------------------------------------------------------------------------------
- def opt_SeqDesign(self, Model, sigma2Dict, NCandidate=5, var='DKL'):
-
- # Initialization
- Model = self.ModelObj
- Bounds = self.BoundTuples
- NrofNewSample = self.ExpDesign.NrofNewSample
- ExploreMethod = self.ExpDesign.ExploreMethod
- ExploitMethod = self.ExpDesign.ExploitMethod
- NrofCandGroups = self.ExpDesign.NrofCandGroups
- TradeOffScheme = self.ExpDesign.TradeOffScheme
-
- OldExpDesign = self.ExpDesign.X
- OutputDictY = self.ExpDesign.Y.copy()
-
- ndim = self.ExpDesign.X.shape[1]
-
- # -----------------------------------------
- # ----------- CUSTOMIZED METHODS ----------
- # -----------------------------------------
- # Utility function ExploitMethod provided by user
- if ExploitMethod.lower() == 'user':
-
- Xnew, filteredSamples = self.ExpDesign.ExploitFunction(self)
-
- print("\n")
- print("\nXnew:\n", Xnew)
-
- return Xnew, filteredSamples
- # -----------------------------------------
- # ---------- EXPLORATION METHODS ----------
- # -----------------------------------------
-
- if ExploreMethod == 'dual annealing':
- # ------- EXPLORATION: OPTIMIZATION -------
- import time
- start_time = time.time()
-
- #Divide the domain to subdomains
- Subdomains = self.subdomain(Bounds, NrofNewSample)
-
- args = [(ExploitMethod, Subdomains[i], sigma2Dict, var, i) for i in range(NrofNewSample)]
-
- # Multiprocessing
- pool = multiprocessing.Pool(multiprocessing.cpu_count())
-
- # With Pool.starmap_async()
- results = pool.starmap_async(self.dual_annealing, args).get()
-
- # Close the pool
- pool.close()
-
- Xnew = np.array([results[NofE][1] for NofE in range(NrofNewSample)])
-
- print("\nXnew:\n", Xnew)
-
- elapsed_time = time.time() - start_time
- print("\n")
- print("elapsed_time: %s sec."%round(elapsed_time,2))
- print('-'*20)
-
- elif ExploreMethod == 'LOOCV':
- # -----------------------------------------------------------------
- # TODO: LOOCV model construnction based on Feng et al. (2020)
- #'LOOCV':
- # Initilize the ExploitScore array
- OutputNames = list(OutputDictY.keys())[1:]
-
- # Generate random samples
- allCandidates = self.ExpDesign.GetSample(NCandidate,
- SamplingMethod='random')
-
- # Construct error model based on LCerror
- errorModel = self.create_ModelError(OldExpDesign, self.LCerror)
- self.errorModel.append(copy.copy(errorModel))
-
- # Evaluate the error models for allCandidates
- eLCAllCands, _ = errorModel.eval_errormodel(allCandidates)
-
- # Select the maximum as the representative error
- eLCAllCandidates = np.max(np.dstack(eLCAllCands.values()), axis=1)[:,0]
-
- # Normalize the error w.r.t the maximum error
- scoreExploration = eLCAllCandidates / np.sum(eLCAllCandidates)
-
- else:
- # ------- EXPLORATION: SPACE-FILLING DESIGN -------
- # Generate candidate samples from Exploration class
- exploration = Exploration(self,NCandidate)
- exploration.w = 100# * ndim #500
- exploration.mcCriterion = 'mc-intersite-proj' #'mc-intersite-proj-th'
- allCandidates , scoreExploration = exploration.getExplorationSamples()
-
- # Temp: ---- Plot all candidates -----
- if ndim == 2:
- def plotter(points, allCandidates, Method, scoreExploration=None):
-
- if Method == 'Voronoi':
- from scipy.spatial import Voronoi, voronoi_plot_2d
- vor = Voronoi(points)
- fig = voronoi_plot_2d(vor)
- ax1 = fig.axes[0]
- else:
- fig = plt.figure()
- ax1 = fig.add_subplot(111)
- ax1.scatter(points[:,0], points[:,1], s=10, c='r', marker="s", label='Old Design Points')
- ax1.scatter(allCandidates[:,0],allCandidates[:,1], s=10, c='b', marker="o", label='Design candidates')
- for i in range(points.shape[0]):
- txt = 'p'+str(i+1)
- ax1.annotate(txt, (points[i,0],points[i,1]))
- if scoreExploration is not None:
- for i in range(allCandidates.shape[0]):
- txt = str(round(scoreExploration[i],5))
- ax1.annotate(txt, (allCandidates[i,0],allCandidates[i,1]))
-
- plt.xlim(self.BoundTuples[0])
- plt.ylim(self.BoundTuples[1])
- # plt.show()
- plt.legend(loc='upper left');
- if ndim == 3:
- def plotter3D(points, closestPoints, Method, scoreExploration=None):
- from mpl_toolkits.mplot3d import Axes3D
-
- jet= plt.get_cmap('Dark2')
- colors = iter(jet(np.linspace(0,1,10)))
-
- fig = plt.figure()
- ax1 = Axes3D(fig)
-
- for idx in range(len(closestPoints)):
- Color = next(colors)
- ax1.scatter(points[idx,0], points[idx,1], points[idx,2], s=10, color = Color, marker="s")
- ax1.scatter(closestPoints[idx][:,0],closestPoints[idx][:,1], closestPoints[idx][:,2], s=10, color = Color, marker="o")
-
- for i in range(points.shape[0]):
- label = 'p'+str(i+1)
- ax1.text(points[i,0],points[i,1],points[i,2], label, None)
-
- if scoreExploration is not None:
- for i in range(allCandidates.shape[0]):
- label = str(round(scoreExploration[i],5))
- ax1.text(allCandidates[i,0],allCandidates[i,1],allCandidates[i,2],label, None)
-
- #ax1.set_xlim(-5, 5)
- #ax1.set_ylim(-5, 5)
- #ax1.set_zlim(-5, 5)
- plt.show()
- plt.legend(loc='upper left');
-
- #plotter3D(self.ExpDesign.X, exploration.closestPoints, ExploreMethod)
-
-
- # -----------------------------------------
- # --------- EXPLOITATION METHODS ----------
- # -----------------------------------------
-
- if ExploitMethod == 'BayesOptDesign' or ExploitMethod == 'BayesActDesign':
-
- # ------- Calculate Exoploration weight -------
- # Compute exploration weight based on trade off scheme
- weightExploration = self.tradoffWeights(TradeOffScheme, OldExpDesign, OutputDictY)
- print("\nweightExploration={0:0.3f} weightExploitation={1:0.3f}".format(weightExploration, (1 - weightExploration)))
-
- # ------- EXPLOITATION: BayesOptDesign & ActiveLearning -------
- if weightExploration != 1.0:
-
- # Create a sample pool for rejection sampling
- MCsize = 15000
- X_MC = self.ExpDesign.GetSample(MCsize, 'random');
-
- # Multiprocessing
- pool = multiprocessing.Pool(multiprocessing.cpu_count())
-
- # Split the candidates in groups for multiprocessing
- split_allCandidates = np.array_split(allCandidates, NrofCandGroups, axis=0)
- args = [(ExploitMethod, Model, split_allCandidates[i], i, sigma2Dict , var, X_MC) for i in range(NrofCandGroups)]
-
- # With Pool.starmap_async()
- results = pool.starmap_async(self.Utility_runner, args).get()
-
- # Close the pool
- pool.close()
-
- # Retrieve the results and append them
- U_J_d = np.concatenate([results[NofE][1] for NofE in range(NrofCandGroups)])
-
- # TODO: Get the expected value (mean) of the Utility score for each cell
- if ExploreMethod == 'Voronoi': U_J_d = np.mean(U_J_d.reshape(-1, NCandidate), axis=1)
-
- # 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)))
- # ------- Calculate Total score -------
- # if ExploreMethod == 'Voronoi':
- # # Randomly select points in each vornoi
- # badSampleIdx=[]
- # nsurvivedOldDesign = len(scoreExploration)
- # allCandidates = np.empty((0, self.ExpDesign.X.shape[1]))
- # for idx in range(nsurvivedOldDesign):
- # candidates = exploration.closestPoints[idx]
- # if len(candidates) == 0:
- # badSampleIdx.append(idx)
- # else:
- # randomIdx = np.random.randint(0, len(candidates), NCandidate)
- # allCandidates = np.vstack((allCandidates, candidates[randomIdx]))
-
- # scoreExploration = np.delete(scoreExploration, badSampleIdx, axis=0)
-
- # ------- Calculate Total score -------
- # ------- Trade off between EXPLORATION & EXPLOITATION -------
- # Total score
- totalScore = (1 - weightExploration)*norm_U_J_d + weightExploration*scoreExploration
-
- # temp: Plot
- # dim = self.ExpDesign.X.shape[1]
- # if dim == 2:
- # plotter(self.ExpDesign.X, allCandidates, ExploreMethod)
-
- # ------- Select the best candidate -------
- # find an optimal point subset to add to the initial design
- # by maximization of the utility score and taking care of NaN values
- temp = totalScore.copy()
- temp[np.isnan(totalScore)] = -np.inf
- sorted_idxtotalScore = np.argsort(temp)[::-1]
- bestIdx = sorted_idxtotalScore[:NrofNewSample]
-
- # select the requested number of samples
- if ExploreMethod == 'Voronoi':
- Xnew = np.zeros((NrofNewSample,ndim))
- for i, idx in enumerate(bestIdx):
- X_can = exploration.closestPoints[idx]
-
- # Calculate the maxmin score for the region of interest
- newSamples, maxminScore = exploration.getMCSamples(allCandidates=X_can)
-
- # select the requested number of samples
- Xnew[i] = newSamples[np.argmax(maxminScore)]
- else:
- Xnew = allCandidates[sorted_idxtotalScore[:NrofNewSample]]
-
- elif ExploitMethod == 'VarOptDesign':
- # ------- EXPLOITATION: VarOptDesign -------
- UtilMethod = var
-
- # ------- Calculate Exoploration weight -------
- # Compute exploration weight based on trade off scheme
- weightExploration = self.tradoffWeights(TradeOffScheme, OldExpDesign, OutputDictY)
- print("\nweightExploration={0:0.3f} weightExploitation={1:0.3f}".format(weightExploration, (1 - weightExploration)))
-
- # Generate candidate samples from Exploration class
- OutputNames = list(OutputDictY.keys())[1:]
- nMeasurement = OutputDictY[OutputNames[0]].shape[1]
-
- # Find indices of the Vornoi cells with samples
- goodSampleIdx=[]
- for idx in range(len(exploration.closestPoints)):
- if len(exploration.closestPoints[idx]) != 0: goodSampleIdx.append(idx)
-
-
- # Find sensitive region
- if UtilMethod == 'LOOCV':
-
- allModifiedLOO = np.zeros((len(OldExpDesign),len(OutputNames), nMeasurement))
- for y_idx , y_key in enumerate(OutputNames):
- for idx, key in enumerate(self.LCerror[y_key].keys()):
- allModifiedLOO[:, y_idx,idx] = abs(self.LCerror[y_key][key])
-
- ExploitScore = np.max(np.max(allModifiedLOO, axis=1),axis=1)
-
- elif UtilMethod in ['EIGF', 'Entropy']:
- # ----- All other in ['EIGF', 'Entropy', 'ALM'] -----
- # Initilize the ExploitScore array
- ExploitScore = np.zeros((len(OldExpDesign),len(OutputNames)))
-
- # Split the candidates in groups for multiprocessing
- if ExploreMethod != 'Voronoi':
- split_allCandidates = np.array_split(allCandidates, NrofCandGroups, axis=0)
- goodSampleIdx = range(NrofCandGroups)
- else:
- split_allCandidates = exploration.closestPoints
-
-
- # Split the candidates in groups for multiprocessing
- args = [(ExploitMethod, Model, split_allCandidates[index], index, sigma2Dict , var) for index in goodSampleIdx]
- #args = [(ExploitMethod, Model, exploration.closestPoints[index], index, sigma2Dict , var) for index in goodSampleIdx]
-
- # Multiprocessing
- pool = multiprocessing.Pool(multiprocessing.cpu_count())
-
- # With Pool.starmap_async()
- results = pool.starmap_async(self.Utility_runner, args).get()
-
- # Close the pool
- pool.close()
-
- # Retrieve the results and append them
- ExploitScore = np.concatenate([results[NofE][1] for NofE in range(len(goodSampleIdx))])
-
- else:
-
- raise NameError('The utility function you requested is not available.')
-
- #print("ExploitScore:\n", ExploitScore)
-
- # find an optimal point subset to add to the initial design
- # by maximization of the utility score and taking care of NaN values
- # Total score
- # Normalize U_J_d
- ExploitScore = ExploitScore / np.sum(ExploitScore)
- totalScore = (1 - weightExploration)*ExploitScore + weightExploration*scoreExploration
-
- temp = totalScore.copy()
- sorted_idxtotalScore = np.argsort(temp, axis=0)[::-1]
- bestIdx = sorted_idxtotalScore[:NrofNewSample]
-
- Xnew = np.zeros((NrofNewSample,ndim))
- if ExploreMethod != 'Voronoi':
- Xnew = allCandidates[bestIdx]
- else:
- for i, idx in enumerate(bestIdx.flatten()):
- X_can = exploration.closestPoints[idx]
-
- #plotter(self.ExpDesign.X, X_can, ExploreMethod, scoreExploration=None)
-
- # Calculate the maxmin score for the region of interest
- newSamples, maxminScore = exploration.getMCSamples(allCandidates=X_can)
-
- # select the requested number of samples
- Xnew[i] = newSamples[np.argmax(maxminScore)]
-
- elif ExploitMethod == 'alphabetic':
- # ------- EXPLOITATION: ALPHABETIC -------
- Xnew = self.util_AlphOptDesign(Model, allCandidates, var)
-
- elif ExploitMethod == 'Space-filling':
- # ------- EXPLOITATION: SPACE-FILLING -------
- totalScore = scoreExploration
-
- # ------- Select the best candidate -------
- # find an optimal point subset to add to the initial design
- # by maximization of the utility score and taking care of NaN values
- temp = totalScore.copy()
- temp[np.isnan(totalScore)] = -np.inf
- sorted_idxtotalScore = np.argsort(temp)[::-1]
-
- # select the requested number of samples
- Xnew = allCandidates[sorted_idxtotalScore[:NrofNewSample]]
-
- else:
- raise NameError('The method you requested is not available.')
-
- print("\n")
- print("\nRun No. {}:".format(OldExpDesign.shape[0]+1))
- print("Xnew:\n", Xnew)
-
- return Xnew, None
-
- #--------------------------------------------------------------------------------------------------------
- def util_AlphOptDesign(self, Model, allCandidates, var='D-Opt'):
- """
- Enrich the Experimental design with the requested alphabetic criterion
- based on exploring the space with number of sampling points (candidates).
-
- Ref: Hadigol, M., & Doostan, A. (2018). Least squares polynomial chaos expansion:
- A review of sampling strategies., Computer Methods in Applied Mechanics and Engineering, 332, 382-407.
-
- Arguments
- ---------
- Model : obj
- Model specification
-
- NCandidate : int
- Number of candidate points to be searched
-
- var : string
- Alphabetic optimality criterion
-
- Returns
- -------
- X_new : ndarray[1, n]
- The new sampling location in the input space.
- """
- NrofNewSample = self.ExpDesign.NrofNewSample
- index = self.index
- NCandidate = allCandidates.shape[0]
-
- # TODO: Loop over outputs
- OutputName = Model.Output.Names[0]
-
- # To avoid changes ub original aPCE object
- from copy import deepcopy
- PCEModel = deepcopy(self)
-
- # Old Experimental design
- oldExpDesignX = PCEModel.ExpDesign.X
-
- # TODO: Only one psi can be selected.
- # Suggestion: Go for the one with the highest LOO error
- Scores = list(PCEModel.ScoresDict[OutputName].values())
- ModifiedLOO = [1-score for score in Scores]
- outIdx = np.argmax(ModifiedLOO[:index])
-
-
- # Initialize Phi to save the criterion's values
- Phi = np.zeros((NCandidate))
-
- BasisIndices = self.BasisDict[OutputName]["y_"+str(outIdx+1)]
- P = len(BasisIndices)
-
-
- # ------ Old Psi ------------
- univ_p_val = self.univ_basis_vals(oldExpDesignX)
- Psi = self.PCE_create_Psi(BasisIndices,univ_p_val)
-
-
- # ------ New candidates (Psi_c) ------------
- # Assemble Psi_c
- univ_p_val_c = self.univ_basis_vals(allCandidates)
- Psi_c = self.PCE_create_Psi(BasisIndices,univ_p_val_c)
-
-
- for idx in range(NCandidate):
-
- # Include the new row to the original Psi
- Psi_cand = np.vstack((Psi, Psi_c[idx]))
-
- # Information matrix
- PsiTPsi = np.dot(Psi_cand.T, Psi_cand)
- M = PsiTPsi / (len(oldExpDesignX)+1)
-
- if np.linalg.cond(PsiTPsi) > 1e-12 and np.linalg.cond(PsiTPsi) < 1/sys.float_info.epsilon:
- # faster
- invM = sp.linalg.solve(M, sp.sparse.eye(PsiTPsi.shape[0]).toarray())
- else:
- # stabler
- invM = np.linalg.pinv(M)
-
- # ---------- Calculate optimality criterion ----------
- # Optimality criteria according to Section 4.5.1 in Ref.
-
- # D-Opt
- if var == 'D-Opt':
- Phi[idx] = (np.linalg.det(invM)) ** (1/P)
-
- # A-Opt
- elif var == 'A-Opt':
- Phi[idx] = np.trace(invM)
-
- # K-Opt
- elif var == 'K-Opt':
- Phi[idx] = np.linalg.cond(M)
-
- else:
- print('The optimality criterion you requested has not been implemented yet!')
-
-
- # find an optimal point subset to add to the initial design
- # by minimization of the Phi
- sorted_idxtotalScore = np.argsort(Phi)
-
- # select the requested number of samples
- Xnew = allCandidates[sorted_idxtotalScore[:NrofNewSample]]
-
- return Xnew
-
- # -------------------------------------------------------------------------
- def eval_metamodel(self,samples=None, nsamples=None, Y=None, discrepModel=False,samplingMethod='random', name='Calib', return_samples=False):
- ModelDict = self.gp_poly if self.metaModel.lower() == 'gpe' else self.CoeffsDict
+ return meanPCEOutputs, stdPCEOutputs
+
+ # -------------------------------------------------------------------------
+ def eval_metamodel(self, samples=None, nsamples=None, Y=None,
+ discrepModel=False, samplingMethod='random',
+ name='Calib', return_samples=False):
+ if self.metaModel.lower() == 'gpe':
+ ModelDict = self.gp_poly
+ else:
+ ModelDict = self.CoeffsDict
if samples is None:
if nsamples is None:
@@ -1928,7 +1045,7 @@ class Metamodel:
else:
self.NrofSamples = nsamples
- samples = self.ExpDesign.GetSample(self.NrofSamples,
+ samples = self.ExpDesign.GetSample(self.NrofSamples,
samplingMethod,
transform=True)
else:
@@ -1976,16 +1093,6 @@ class Metamodel:
y_mean = np.dot(psi, coeffs)
y_std = np.zeros_like(y_mean)
- # if self.metaModel.lower() == 'pcekriging':
- # gp = self.gp_poly[Outkey][Inkey]
- # y_gp_mean, y_gp_std = gp.predict(samples, return_std=True)
- # y_mean += y_gp_mean
- # try:
- # y_std = y_error[Outkey][:,Inkey]
- # print("y_std:",y_std)
- # except:
- # y_std = np.zeros(shape=y_mean.shape)
-
PCEOutputs_mean[:, idx] = y_mean
PCEOutputs_std[:, idx] = y_std
idx += 1
@@ -1993,8 +1100,11 @@ class Metamodel:
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))
+ meanPCEOutputs[Outkey] = PCA.mean_
+ meanPCEOutputs[Outkey] += 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
@@ -2004,16 +1114,22 @@ class Metamodel:
if 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]
+ meanPCEOutputs[Outkey] = PCA.mean_[:index]
+ meanPCEOutputs[Outkey] += 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:]
+ meanPCEOutputs[Outkey] = PCA.mean_[index:]
+ meanPCEOutputs[Outkey] += 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:]
@@ -2023,1012 +1139,3 @@ class Metamodel:
return meanPCEOutputs, stdPCEOutputs, samples
else:
return meanPCEOutputs, stdPCEOutputs
-
- #--------------------------------------------------------------------------------------------------------
- def normpdf(self, PCEOutputs, std_PC_MC, ObservationData, Sigma2s):
-
- ModelOutputNames = self.ModelObj.Output.Names
-
- SampleSize, index = PCEOutputs[ModelOutputNames[0]].shape
- # if self.index is not None: index = self.index
- index = self.index if type(self.index) is list else [self.index]*len(ModelOutputNames)
-
- # Flatten the ObservationData
- TotalData = ObservationData[ModelOutputNames].to_numpy().flatten('F')
-
- # Remove NaN
- TotalData = TotalData[~np.isnan(TotalData)]
- Sigma2s = Sigma2s[~np.isnan(Sigma2s)]
-
- # Flatten the Output
- TotalOutputs = np.empty((SampleSize,0))
- for idx, key in enumerate(ModelOutputNames):
- TotalOutputs = np.hstack((TotalOutputs, PCEOutputs[key][:,:index[idx]]))
-
- # Covariance Matrix
- covMatrix = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
- np.fill_diagonal(covMatrix, Sigma2s)
-
- # Add the std of the PCE.
- covMatrix_PCE = np.zeros((Sigma2s.shape[0], Sigma2s.shape[0]), float)
- stdPCE = np.empty((SampleSize,0))
- for idx, key in enumerate(ModelOutputNames):
- stdPCE = np.hstack((stdPCE, std_PC_MC[key][:,:index[idx]]))
-
- # Expected value of variance (Assump: i.i.d stds)
- varPCE = np.mean(stdPCE**2, axis=0)
- # # varPCE = np.var(stdPCE, axis=1)
- np.fill_diagonal(covMatrix_PCE, varPCE)
-
- # Aggregate the cov matrices
- covMatrix += covMatrix_PCE
-
- # Compute likelihood
- self.Likelihoods = stats.multivariate_normal.pdf(TotalOutputs,
- mean=TotalData,
- cov=covMatrix,
- allow_singular=True)
-
- return self.Likelihoods
-
- #--------------------------------------------------------------------------------------------------------
- def BME_Corr_Weight(self, PCEModel, Data, Sigma2sDict):
- """
- Calculates the correction factor for BMEs.
- """
-
- # ----- Original model outputs on the ExpDesign ---------------
- # Extract the outputs
- Model = PCEModel.ModelObj
- OutputType = Model.Output.Names
- OrigModelOutput = PCEModel.ExpDesign.Y
- SampleSize = OrigModelOutput[OutputType[0]].shape[0]
-
- # Flatten the OutputType for OrigModel
- TotalOutputs = np.empty((SampleSize,0))
- for outputType in OutputType:
- TotalOutputs = np.hstack((TotalOutputs, OrigModelOutput[outputType]))
-
- # ----- PCE model outputs on the ExpDesign ---------------
- # Evaluate PCEModel on the experimental design
- Samples = PCEModel.ExpDesign.X
- OutputRS, stdOutputRS = PCEModel.eval_metamodel(samples=Samples)
-
- # Flatten the mean Outputs for PCEModel
- TotalPCEOutputs = np.empty((SampleSize,0))
- for outputType in OutputType:
- TotalPCEOutputs = np.hstack((TotalPCEOutputs, OutputRS[outputType]))
-
- NofMeasurements = TotalPCEOutputs.shape[1]
-
- # Flatten the mean Outputs for PCEModel
- stdPCEOutputs = np.empty((SampleSize,0))
- for outputType in OutputType:
- stdPCEOutputs = np.hstack((stdPCEOutputs, stdOutputRS[outputType]))
-
- # ----- Observation data ---------------
- # Flatten Oberservation datasets
- #Observations = Data.flatten()
- Observations = np.empty((1,0))
- for outputType in OutputType:
- Observations = np.append(Observations, Data[outputType])
-
- # ------------- Sigma2s ---------------
- Sigma2s = np.empty((1,0))
- for outputType in OutputType:
- Sigma2s = np.append(Sigma2s, Sigma2sDict[outputType])
-
- # ----------------- LIKELIHOOD CALCULATION --------------------
- # Initialization
- covMatrix=np.zeros((NofMeasurements, NofMeasurements), float)
- BME_RM_Model_Weight = np.zeros((SampleSize))
- BME_RM_Data_Weight = np.zeros((SampleSize))
- BME_Corr = 0
-
- # Deviation Computations
- RM_Model_Deviation = np.zeros((SampleSize,NofMeasurements))
- RM_Data_Deviation = np.zeros((SampleSize,NofMeasurements))
- for i in range(SampleSize):
- RM_Model_Deviation[i] = TotalOutputs[i][:NofMeasurements] - TotalPCEOutputs[i, :] # Reduce model- Full Model
- RM_Data_Deviation[i] = Observations - TotalPCEOutputs[i, :] # Reduce model- Measurement Data
-
-
- # Initialization of Co-Variance Matrix
- # For BME_RM_ModelWeight
- if NofMeasurements == 1:
- RM_Model_Error = np.zeros((NofMeasurements, NofMeasurements), float)
- np.fill_diagonal(RM_Model_Error, np.cov(RM_Model_Deviation.T))
- else:
- RM_Model_Error = np.cov(RM_Model_Deviation.T)
-
-
- # Computation of Weight according to the deviations
- for i in range(SampleSize):
- # For BME_RM_DataWeight
- try:
- var = Sigma2s[i]
- if len(var)==1:
- np.fill_diagonal(covMatrix, var)
- else:
- row,col = np.diag_indices(covMatrix.shape[0])
- covMatrix[row,col] = np.hstack((np.repeat(var[0], NofMeasurements*0.5),np.repeat(var[1], NofMeasurements*0.5)))
-
- except:
- var = Sigma2s
-
- np.fill_diagonal(covMatrix, var)
-
- # Add the std of the PCE is emulator is chosen.
- covMatrix_PCE = np.zeros((NofMeasurements, NofMeasurements), float)
-
- stdPCE = np.mean(stdPCEOutputs, axis=1)
- np.fill_diagonal(covMatrix_PCE, stdPCE**2)
-
- covMatrix = covMatrix + covMatrix_PCE
-
- # Calculate the denomitor
- denom1 = (np.sqrt(2*np.pi)) ** NofMeasurements
- denom2 = (((2*np.pi)**(NofMeasurements/2)) * np.sqrt(np.linalg.det(covMatrix)))
-
- BME_RM_Model_Weight[i] = (np.exp(-0.5 * np.dot(np.dot(RM_Model_Deviation[i], np.linalg.pinv(RM_Model_Error)), RM_Model_Deviation[i])))/denom1
- BME_RM_Data_Weight[i] = (np.exp(-0.5 * np.dot(np.dot(RM_Data_Deviation[i], np.linalg.pinv(covMatrix)), RM_Data_Deviation[i][:,np.newaxis])))/denom2
-
- for i in range(SampleSize):
- BME_Corr += BME_RM_Model_Weight[i] * BME_RM_Data_Weight[i] / np.nansum(BME_RM_Data_Weight)
-
-
- return np.log(BME_Corr)
-
- #--------------------------------------------------------------------------------------------------------
- def posteriorPlot(self, Posterior, MAP, parNames, key, figsize=(10,10)):
-
- # Initialization
- newpath = (r'Outputs_SeqPosteriorComparison')
- if not os.path.exists(newpath): os.makedirs(newpath)
-
- lw = 3.
- BoundTuples = self.BoundTuples
- NofPa = len(MAP)
-
- # This is the true mean of the second mode that we used above:
- value1 = MAP
-
- # This is the empirical mean of the sample:
- value2 = np.mean(Posterior, axis=0)
-
- if NofPa == 2:
-
- figPosterior, ax = plt.subplots()
- plt.hist2d(Posterior[:,0], Posterior[:,1], bins=(200, 200),
- range=np.array([BoundTuples[0],BoundTuples[1]]), cmap=plt.cm.jet)
-
- plt.xlabel(parNames[0])
- plt.ylabel(parNames[1])
-
- plt.xlim(BoundTuples[0])
- plt.ylim(BoundTuples[1])
-
- ax.axvline(value1[0], color="g", lw=lw)
- ax.axhline(value1[1], color="g", lw=lw)
- ax.plot(value1[0], value1[1], "sg", lw=lw+1)
-
- ax.axvline(value2[0], ls='--', color="r", lw=lw)
- ax.axhline(value2[1], ls='--', color="r", lw=lw)
- ax.plot(value2[0], value2[1], "sr", lw=lw+1)
-
- else:
- import corner
- figPosterior = corner.corner(Posterior, labels=parNames,title_fmt='.2e',
- show_titles=True, title_kwargs={"fontsize": 12})
-
- # Extract the axes
- axes = np.array(figPosterior.axes).reshape((NofPa, NofPa))
-
- # Loop over the diagonal
- for i in range(NofPa):
- ax = axes[i, i]
- ax.axvline(value1[i], color="g")
- ax.axvline(value2[i], ls='--', color="r")
-
- # Loop over the histograms
- for yi in range(NofPa):
- for xi in range(yi):
- ax = axes[yi, xi]
- ax.axvline(value1[xi], color="g")
- ax.axvline(value2[xi], ls='--', color="r")
- ax.axhline(value1[yi], color="g")
- ax.axhline(value2[yi], ls='--', color="r")
- ax.plot(value1[xi], value1[yi], "sg")
- ax.plot(value2[xi], value2[yi], "sr")
-
-
- #figPosterior.set_size_inches((10,10))
- figPosterior.savefig('./'+newpath+'/'+key+'.svg',bbox_inches='tight')
- plt.close()
-
- # Save the posterior as .npy
- np.save('./'+newpath+'/'+key+'.npy', Posterior)
-
- return figPosterior
-
- #--------------------------------------------------------------------------------------------------------
- def hellinger_distance(self, P, Q):
- """
- Hellinger distance between two continuous distributions.
-
- The maximum distance 1 is achieved when P assigns probability zero to
- every set to which Q assigns a positive probability, and vice versa.
- 0 (identical) and 1 (maximally different)
- """
-
- mu1 = P.mean()
- Sigma1 = np.std(P)
-
- mu2 = Q.mean()
- Sigma2 = np.std(Q)
-
- term1 = np.sqrt(2*Sigma1*Sigma2 / (Sigma1**2 + Sigma2**2))
-
- term2 = np.exp(-.25 * (mu1 - mu2)**2 / (Sigma1**2 + Sigma2**2))
-
- H_squared = 1 - term1 * term2
-
-
- return np.sqrt(H_squared)
-
- #--------------------------------------------------------------------------------------------------------
- def BME_Calculator(self, PCEModel, ObservationData, sigma2Dict):
- """
- This function computes the Bayesian model evidence (BME) via Monte Carlo
- integration.
-
- """
-
- # Initializations
- validLikelihoods = self.validLikelihoods
- PostSnapshot = self.ExpDesign.PostSnapshot
- if PostSnapshot or len(validLikelihoods)!=0:
- newpath = (r'Outputs_SeqPosteriorComparison')
- if not os.path.exists(newpath): os.makedirs(newpath)
-
- SamplingMethod = 'random'
- MCsize = 100000
- ESS = 0
-
- # Estimation of the integral via Monte Varlo integration
- while ((ESS > MCsize) or (ESS < 1)):
-
- # Generate samples for Monte Carlo simulation
- if len(validLikelihoods)==0:
- # ExpDesign = ExpDesigns(self.Inputs)
- X_MC = PCEModel.ExpDesign.GetSample(MCsize, SamplingMethod);
- else:
- X_MC = self.validSamples
- MCsize = X_MC.shape[0]
-
- # Monte Carlo simulation for the candidate design
- 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)
- Likelihoods = self.normpdf(Y_PC_MC,std_PC_MC,ObservationData, sigma2Dict)
-
- # Check the Effective Sample Size (1000<ESS<MCsize)
- ESS = 1 / np.sum(np.square(Likelihoods/np.sum(Likelihoods)))
-
- # Enlarge sample size if it doesn't fulfill the criteria
- if ((ESS > MCsize) or (ESS < 1)):
- print('ESS={0} MC size should be larger.'.format(ESS))
- MCsize = MCsize * 10
- ESS = 0
-
- # Rejection Step
- # Random numbers between 0 and 1
- unif = np.random.rand(1, MCsize)[0]
-
- # Reject the poorly performed prior
- accepted = (Likelihoods/np.max(Likelihoods)) >= unif
- X_Posterior = X_MC[accepted]
-
-
- # -------------------------------------------------------------------------
- # ------------- Kullback-Leibler Divergence & Information Entropy ---------
- # -------------------------------------------------------------------------
- # Prior-based estimation of BME
- logBME = np.log(np.nanmean(Likelihoods))
- # Calculate BME and its correction factor
-# BMECorrFactor = self.BME_Corr_Weight(PCEModel, ObservationData, sigma2Dict)
-#
-# logBME += BMECorrFactor
-
- # Posterior-based expectation of likelihoods
- postExpLikelihoods = np.mean(np.log(Likelihoods[accepted]))
-
- # Posterior-based expectation of prior densities
- #postExpPrior = np.mean([log_prior(sample) for sample in X_Posterior])
-
-
- # Calculate Kullback-Leibler Divergence
- #KLD = np.mean(np.log(Likelihoods[Likelihoods!=0])- logBME)
- KLD = postExpLikelihoods - logBME
-
- # Information Entropy based on Entropy paper Eq. 38
- #infEntropy = logBME - postExpPrior - postExpLikelihoods
-
- # If PostSnapshot is True, plot likelihood vs refrence
- if PostSnapshot or len(validLikelihoods)!=0:
- idx = len([name for name in os.listdir(newpath) if 'Likelihoods_' in name and os.path.isfile(os.path.join(newpath, name))])
- fig, ax = plt.subplots()
- sns.kdeplot(np.log(validLikelihoods[validLikelihoods>0]), shade=True, color="g", label='Ref. Likelihood')
- sns.kdeplot(np.log(Likelihoods[Likelihoods>0]), shade=True, color="b", label= 'Likelihood with PCE')
-
- # Hellinger distance
- distHellinger = self.hellinger_distance(np.log(validLikelihoods[validLikelihoods>0]), np.log(Likelihoods[Likelihoods>0]))
- text = "Hellinger Dist.={:.3f}\n logBME={:.3f}\n DKL={:.3f}".format(distHellinger, logBME, KLD)
-
- plt.text(0.05, 0.75, text,bbox=dict(facecolor='wheat', edgecolor='black', boxstyle='round,pad=1'),
- transform = ax.transAxes)
- #plt.show()
-
- fig.savefig('./'+newpath+'/Likelihoods_'+str(idx)+'.svg',bbox_inches='tight')
- plt.close()
-
- else:
- distHellinger = 0.0
-
- return (logBME, KLD, X_Posterior,distHellinger)
-
- #--------------------------------------------------------------------------------------------------------
-
- def validError_(self):
-
- Model = self.ModelObj
- OutputName = Model.Output.Names
-
- # Generate random samples
- Samples = self.validSamples
-
- # Extract the original model with the generated samples
- ModelOutputs = self.validModelRuns
-
- # Run the PCE model with the generated samples
- PCEOutputs, PCEOutputs_std = self.eval_metamodel(samples=Samples)
-
- validError_dict = {}
- # Loop over the keys and compute RMSE error.
- for key in OutputName:
- weight = np.mean(np.square(PCEOutputs_std[key]), axis=0)
- if all(weight==0): weight='variance_weighted'
- validError_dict[key] = mean_squared_error(ModelOutputs[key], PCEOutputs[key])/\
- np.var(ModelOutputs[key],ddof=1)
-
- return validError_dict
-
- #--------------------------------------------------------------------------------------------------------
- def error_Mean_Std(self):
-
- from sklearn.metrics import mean_squared_error, mean_absolute_error
- # Extract the mean and std provided by user
- df_MCReference = self.ModelObj.read_MCReference()
-
- # Compute the mean and std based on the PCEModel
- PCEMeans = dict()
- PCEStds = dict()
- for Outkey, ValuesDict in self.CoeffsDict.items():
-
- PCEMean = np.zeros((len(ValuesDict)))
- PCEStd = np.zeros((len(ValuesDict)))
-
- for Inkey, InIdxValues in ValuesDict.items():
- idx = int(Inkey.split('_')[1]) - 1
- coeffs = self.CoeffsDict[Outkey][Inkey]
-
- # Mean = c_0
- PCEMean[idx] = coeffs[0] if coeffs[0]!=0 else self.clf_poly[Outkey][Inkey].intercept_
- # Std = sqrt(sum(coeffs[1:]**2))
- PCEStd[idx] = np.sqrt(np.sum(np.square(coeffs[1:])))
-
- if self.DimRedMethod.lower() == 'pca':
- 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)
- # RMSE_std = mean_absolute_error(df_MCReference['std'], PCEStd)
- RMSE_Mean = mean_squared_error(df_MCReference['mean'], PCEMean,squared=False)
- RMSE_std = mean_squared_error(df_MCReference['std'], PCEStd,squared=False)
-
- PCEMeans[Outkey] = PCEMean
- PCEStds[Outkey] = PCEStd
-
-
- return RMSE_Mean, RMSE_std
- #--------------------------------------------------------------------------------------------------------
- def create_PCE_SeqDesign(self, Model):
-
- #Initialization
- errorIncreases = False
-
- # Get the parameters
- MaxNSamples = self.ExpDesign.MaxNSamples
- ModifiedLOOThreshold = self.ExpDesign.ModifiedLOOThreshold
- NCandidate = self.ExpDesign.NCandidate
- UtilityFunctions = self.ExpDesign.UtilityFunction
- PostSnapshot = self.ExpDesign.PostSnapshot
- nReprications = self.ExpDesign.nReprications
-
- # Handle if only one UtilityFunctions is provided
- if not isinstance(self.ExpDesign.UtilityFunction, list):
- UtilityFunctions = [self.ExpDesign.UtilityFunction]
-
- # ---------- Initial PCEModel ----------
- PCEModel = self.create_PCE_normDesign(Model)
- initPCEModel= copy.deepcopy(PCEModel)
-
- # TODO: Loop over outputs
- OutputName = Model.Output.Names
-
- # Estimation of the integral via Monte Varlo integration
- ObservationData = self.Observations
-
- # Check if data is provided
- TotalSigma2 = np.empty((0,1))
- if len(ObservationData) != 0:
- # ---------- Prepare diagonal enteries for co-variance matrix ----------
- for keyIdx, key in enumerate(Model.Output.Names):
- optSigma = 'B'
- sigma2 = np.array(self.Discrepancy.Parameters[key])
- TotalSigma2 = np.append(TotalSigma2,sigma2)
-
- # Calculate the initial BME
- initBME, initKLD, initPosterior, initDistHellinger = self.BME_Calculator(initPCEModel, ObservationData, TotalSigma2)
- print("\nInitial BME:", initBME)
- print("Initial KLD:", initKLD)
-
- # Posterior snapshot (initial)
- if PostSnapshot:
- MAP = self.ExpDesign.MAP
- parNames = self.ExpDesign.parNames
- print('Posterior snapshot (initial) is being plotted...')
- self.posteriorPlot(initPosterior, MAP, parNames, 'SeqPosterior_init')
-
- # Check the convergence of the Mean&Std
- if self.ModelObj.MCReference and self.metaModel.lower()!='gpe':
- initRMSEMean, initRMSEStd = self.error_Mean_Std()
- print("Initial Mean and Std error: %s, %s"%(initRMSEMean, initRMSEStd))
-
- # Read the initial experimental design
- Xinit = initPCEModel.ExpDesign.X
- initTotalNSamples = len(self.ExpDesign.X)
- initYprev = initPCEModel.ModelOutputDict
- initLCerror = initPCEModel.LCerror
-
- # Read the initial ModifiedLOO
- if self.metaModel.lower()!='gpe':
- Scores_all, varExpDesignY =[], []
- for OutName in OutputName:
- Scores_all.append(list(initPCEModel.ScoresDict[OutName].values()))
- if self.DimRedMethod.lower() == 'pca':
- components = self.pca[OutName].transform(initPCEModel.ExpDesign.Y[OutName])
- varExpDesignY.append(np.var(components,axis=0))
- else:
- 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)]
-
- if len(self.validModelRuns) != 0:
- initValidError = initPCEModel.validError_()
- initValidError = list(initValidError.values())
- print("\nInitial ValidError:", initValidError)
-
-
- # Replicate the sequential design
- for repIdx in range(nReprications):
- print ('>>>> Replication: %s<<<<'%(repIdx+1))
-
- # To avoid changes ub original aPCE object
- # PCEModel = copy.deepcopy(initPCEModel)
- PCEModel.ExpDesign.X = Xinit
- PCEModel.ExpDesign.Y = initYprev
- PCEModel.LCerror = initLCerror
-
- for UtilityFunction in UtilityFunctions:
-
- print ('>>>> UtilityFunction: %s<<<<'%UtilityFunction)
- # To avoid changes ub original aPCE object
- # PCEModel = copy.deepcopy(initPCEModel)
- PCEModel.ExpDesign.X = Xinit
- PCEModel.ExpDesign.Y = initYprev
- PCEModel.LCerror = initLCerror
-
-
- # Set the experimental design
- Xprev = Xinit
- TotalNSamples = initTotalNSamples
- Yprev = initYprev
-
- Xfull = []
- Yfull = []
-
- # Store the initial ModifiedLOO
- if self.metaModel.lower()!='gpe':
- print("\nInitial ModifiedLOO:", initModifiedLOO)
-
- ModifiedLOO = initModifiedLOO
- SeqModifiedLOO = np.array(ModifiedLOO)
-
- if len(self.validModelRuns) != 0:
- ValidError = initValidError
- SeqValidError = np.array(ValidError)
-
- # Check if data is provided
- if len(ObservationData) != 0:
- SeqBME = np.array([initBME])
- SeqKLD = np.array([initKLD])
- SeqDistHellinger = np.array([initDistHellinger])
- if self.ModelObj.MCReference and self.metaModel.lower()!='gpe':
- seqRMSEMean = np.array([initRMSEMean])
- seqRMSEStd = np.array([initRMSEStd])
-
- postcnt = 1
- itrNr = 1
- # Start Sequential Experimental Design
- while (TotalNSamples < MaxNSamples) and any(LOO > ModifiedLOOThreshold for LOO in ModifiedLOO):
-
- # Optimal Bayesian Design
- PCEModel.ExpDesignFlag = 'sequential'
- Xnew, updatedPrior = PCEModel.opt_SeqDesign(Model, TotalSigma2, NCandidate, UtilityFunction)
-
-
- S = np.min(distance.cdist(Xinit, Xnew, 'euclidean'))
- self.seqMinDist.append(S)
- print("\nmin Dist from OldExpDesign:", S)
- print("\n")
-
- # Evaluate the full model response at the new sample points:
- Ynew, _ = Model.Run_Model_Parallel(Xnew, prevRun_No=TotalNSamples)
- TotalNSamples += Xnew.shape[0]
-
- # -------- Plot the surrogate model vs Origninal Model -------
- if self.adaptVerbose:
- from PostProcessing.adaptPlot import adaptPlot
- y_hat, std_hat = PCEModel.eval_metamodel(samples=Xnew)
- adaptPlot(PCEModel, Ynew, y_hat, std_hat,plotED=False)
-
- # -------- Retrain the surrogate model -------
- # Extend new experimental design
- Xfull = np.vstack((Xprev,Xnew))
-
- # Updating existing key's value
- for OutIdx in range(len(OutputName)):
- OutName = OutputName[OutIdx]
- try:
- Yfull = np.vstack((Yprev[OutName],Ynew[OutName]))
- except:
- Yfull = np.vstack((Yprev[OutName],Ynew))
-
- PCEModel.ModelOutputDict[OutName] = Yfull
-
-
- PCEModel.ExpDesign.SamplingMethod = 'user'
- PCEModel.ExpDesign.X = Xfull
- PCEModel.ExpDesign.Y = PCEModel.ModelOutputDict
-
- # save the Experimental Design for next iteration
- Xprev = Xfull
- Yprev = PCEModel.ModelOutputDict
-
- # Pass the new prior as the input
- PCEModel.Inputs.polycoeffsFlag = False
- if updatedPrior is not None:
- PCEModel.Inputs.polycoeffsFlag = True
- print("updatedPrior:", updatedPrior.shape)
- # arbitrary polynomial chaos
- for i in range(updatedPrior.shape[1]):
- PCEModel.Inputs.Marginals[i].DistType = None
- PCEModel.Inputs.Marginals[i].InputValues = updatedPrior[:,i]
-
- prevPCEModel = PCEModel
- PCEModel = PCEModel.create_PCE_normDesign(Model)
-
- # -------- Evaluate the retrain surrogate model -------
- # Compute the validation error
- if len(self.validModelRuns) != 0:
- validError = PCEModel.validError_()
- ValidError = list(validError.values())
- print("\nUpdated ValidError:", ValidError)
-
- # Extract Modified LOO from Output
- if self.metaModel.lower()!='gpe':
- Scores_all, varExpDesignY =[], []
- for OutName in OutputName:
- Scores_all.append(list(PCEModel.ScoresDict[OutName].values()))
- if self.DimRedMethod.lower() == 'pca':
- components = self.pca[OutName].transform(initPCEModel.ExpDesign.Y[OutName])
- varExpDesignY.append(np.var(components,axis=0))
- else:
- 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],
- weights=weights)]
-
- print('\n')
- print("Updated ModifiedLOO %s:\n"%UtilityFunction, ModifiedLOO)
- print("Xfull:", Xfull.shape)
- print('\n')
-
-
- # check the direction of the error (on average):
- # if it increases consistently stop the iterations
- n_checks = 3
- if itrNr > n_checks * self.ExpDesign.NrofNewSample:
- ss = np.sign(SeqModifiedLOO - ModifiedLOO) #ss<0 error increasing
- errorIncreases = np.sum(np.mean(ss[-2:], axis=1)) <= -1*n_checks
-
- # If error is increasing in the last n_check iteration,
- # stop the search and return the previous PCEModel
- if errorIncreases:
- print("Warning: The modified error is increasing compared to \
- the last %s iterations."%n_checks)
-
- PCEModel = prevPCEModel
- break
-
- else:
- prevPCEModel = PCEModel
-
-
- # Store updated ModifiedLOO
- SeqModifiedLOO = np.vstack((SeqModifiedLOO, ModifiedLOO))
- if len(self.validModelRuns) != 0:
- SeqValidError = np.vstack((SeqValidError, ValidError))
-
- # -------- Caclulation of BME as accuracy metric -------
- # Check if data is provided
- if len(ObservationData) != 0:
- # Calculate the initial BME
- BME, KLD, Posterior, DistHellinger = self.BME_Calculator(PCEModel, ObservationData, TotalSigma2)
- print('\n')
- print("Updated BME:", BME)
- print("Updated KLD:", KLD)
- print('\n')
-
- # Plot some snapshots of the posterior
- stepSnapshot = self.ExpDesign.stepSnapshot
- if PostSnapshot and postcnt % stepSnapshot == 0:
- MAP = self.ExpDesign.MAP
- parNames = self.ExpDesign.parNames
- print('Posterior snapshot is being plotted...')
- self.posteriorPlot(Posterior, MAP, parNames, 'SeqPosterior_'+str(postcnt))
- postcnt += 1
-
- # Check the convergence of the Mean&Std
- if self.ModelObj.MCReference and self.metaModel.lower()!='gpe':
- print('\n')
- RMSE_Mean, RMSE_std = self.error_Mean_Std()
- print("Updated Mean and Std error: %s, %s"%(RMSE_Mean, RMSE_std))
- print('\n')
-
- # Store the updated BME & KLD &
- # Check if data is provided
- if len(ObservationData) != 0:
- SeqBME = np.vstack((SeqBME, BME))
- SeqKLD = np.vstack((SeqKLD, KLD))
- SeqDistHellinger= np.vstack((SeqDistHellinger, DistHellinger))
- if self.ModelObj.MCReference and self.metaModel.lower()!='gpe':
- seqRMSEMean = np.vstack((seqRMSEMean, RMSE_Mean))
- seqRMSEStd = np.vstack((seqRMSEStd, RMSE_std))
-
- itrNr += 1
- # Store updated ModifiedLOO and BME in dictonary
- strKey = UtilityFunction + '_rep_' +str(repIdx+1)
- PCEModel.SeqModifiedLOO[strKey] = SeqModifiedLOO
- if len(self.validModelRuns) != 0:
- PCEModel.seqValidError[strKey] = SeqValidError
-
- # Check if data is provided
- if len(ObservationData) != 0:
- PCEModel.SeqBME[strKey] = SeqBME
- PCEModel.SeqKLD[strKey] = SeqKLD
- if len(self.validLikelihoods) != 0:
- PCEModel.SeqDistHellinger[strKey] = SeqDistHellinger
- if self.ModelObj.MCReference and self.metaModel.lower()!='gpe':
- PCEModel.seqRMSEMean[strKey] = seqRMSEMean
- PCEModel.seqRMSEStd[strKey] = seqRMSEStd
-
-
- return PCEModel
-
- #--------------------------------------------------------------------------------------------------------
- def create_GPE_SeqDesign(self, Model):
-
- #Initialization
-
- # Get the parameters
- MaxNSamples = self.ExpDesign.MaxNSamples
- NCandidate = self.ExpDesign.NCandidate
- UtilityFunctions = self.ExpDesign.UtilityFunction
- PostSnapshot = self.ExpDesign.PostSnapshot
- nReprications = self.ExpDesign.nReprications
-
- # Handle if only one UtilityFunctions is provided
- if not isinstance(self.ExpDesign.UtilityFunction, list):
- UtilityFunctions = [self.ExpDesign.UtilityFunction]
-
- # ---------- Initial PCEModel ----------
- metamodel = self.create_PCE_normDesign(Model)
- initPCEModel= metamodel
-
- # TODO: Loop over outputs
- OutputName = Model.Output.Names
-
- # Estimation of the integral via Monte Varlo integration
- ObservationData = self.Observations
-
- # Check if data is provided
- TotalSigma2 = np.empty((0,1))
- if len(ObservationData) != 0:
- # ---------- Prepare diagonal enteries for co-variance matrix ----------
- for keyIdx, key in enumerate(Model.Output.Names):
- optSigma = 'B'
- sigma2 = np.array(self.Discrepancy.Parameters[key])
- TotalSigma2 = np.append(TotalSigma2,sigma2)
-
- # Calculate the initial BME
- initBME, initKLD, initPosterior, initDistHellinger = self.BME_Calculator(initPCEModel, ObservationData, TotalSigma2)
- print("\nInitial BME:", initBME)
- print("Initial KLD:", initKLD)
-
- # Posterior snapshot (initial)
- if PostSnapshot:
- MAP = self.ExpDesign.MAP
- parNames = self.ExpDesign.parNames
- print('Posterior snapshot (initial) is being plotted...')
- self.posteriorPlot(initPosterior, MAP, parNames, 'SeqPosterior_init')
-
- # Read the initial experimental design
- Xinit = initPCEModel.ExpDesign.X
- initTotalNSamples = len(self.ExpDesign.X)
- initYprev = initPCEModel.ModelOutputDict
-
- # Read the initial ModifiedLOO
- if len(self.validModelRuns) != 0:
- initValidError = initPCEModel.validError_()
- initValidError = list(initValidError.values())
- print("\nInitial ValidError:", initValidError)
-
- self.SeqModifiedLOO = {}
- self.SeqBME = {}
- self.seqRMSEMean = {}
- self.seqRMSEStd = {}
- self.seqMinDist = []
-
- # Replicate the sequential design
- for repIdx in range(nReprications):
- print ('>>>> Replication: %s<<<<'%(repIdx+1))
- # To avoid changes ub original aPCE object
- metamodel.ExpDesign.X = Xinit
- metamodel.ExpDesign.Y = initYprev
-
- for UtilityFunction in UtilityFunctions:
-
- print ('>>>> UtilityFunction: %s<<<<'%UtilityFunction)
- # To avoid changes ub original aPCE object
- metamodel.ExpDesign.X = Xinit
- metamodel.ExpDesign.Y = initYprev
-
- # Set the experimental design
- Xprev = Xinit
- TotalNSamples = initTotalNSamples
- Yprev = initYprev
-
- Xfull = []
- Yfull = []
-
- # Store the validation error
- if len(self.validModelRuns) != 0:
- ValidError = initValidError
- SeqValidError = np.array(ValidError)
-
- # Check if data is provided
- if len(ObservationData) != 0:
- SeqBME = np.array([initBME])
- SeqKLD = np.array([initKLD])
- SeqDistHellinger = np.array([initDistHellinger])
-
- postcnt = 1
- itrNr = 1
- # Start Sequential Experimental Design
- while TotalNSamples < MaxNSamples:
-
- # Optimal Bayesian Design
- self.ExpDesignFlag = 'sequential'
- Xnew, updatedPrior = self.opt_SeqDesign(Model, TotalSigma2, NCandidate, UtilityFunction)
-
-
- S = np.min(distance.cdist(Xinit, Xnew, 'euclidean'))
- self.seqMinDist.append(S)
- print("\nmin Dist from OldExpDesign:", S)
- print("\n")
-
- # Evaluate the full model response at the new sample points:
- Ynew, _ = Model.Run_Model_Parallel(Xnew, prevRun_No=TotalNSamples)
- TotalNSamples += Xnew.shape[0]
-
- # -------- Plot the surrogate model vs Origninal Model -------
- if self.adaptVerbose:
- from PostProcessing.adaptPlot import adaptPlot
- y_hat, std_hat = metamodel.eval_metamodel(samples=Xnew)
- adaptPlot(metamodel, Ynew, y_hat, std_hat,plotED=True)
-
- # -------- Retrain the surrogate model -------
- # Extend new experimental design
- Xfull = np.vstack((Xprev,Xnew))
-
- # Updating existing key's value
- for OutIdx in range(len(OutputName)):
- OutName = OutputName[OutIdx]
- try:
- Yfull = np.vstack((Yprev[OutName],Ynew[OutName]))
- except:
- Yfull = np.vstack((Yprev[OutName],Ynew))
-
- metamodel.ModelOutputDict[OutName] = Yfull
-
-
- metamodel.ExpDesign.SamplingMethod = 'user'
- metamodel.ExpDesign.X = Xfull
- metamodel.ExpDesign.Y = metamodel.ModelOutputDict
-
- # save the Experimental Design for next iteration
- Xprev = Xfull
- Yprev = metamodel.ModelOutputDict
-
- # Pass the new prior as the input
- if updatedPrior is not None:
- print("updatedPrior:", updatedPrior.shape)
- # arbitrary polynomial chaos
- for i in range(updatedPrior.shape[1]):
- self.Inputs.Marginals[i].DistType = None
- self.Inputs.Marginals[i].InputValues = updatedPrior[:,i]
-
- prevPCEModel = metamodel
- metamodel = metamodel.create_PCE_normDesign(Model)#, OptDesignFlag=True)
-
- # -------- Evaluate the retrain surrogate model -------
- # Compute the validation error
- if len(self.validModelRuns) != 0:
- validError = metamodel.validError_()
- ValidError = list(validError.values())
- print("\nUpdated ValidError:", ValidError)
-
- # Store updated validation error
- if len(self.validModelRuns) != 0:
- SeqValidError = np.vstack((SeqValidError, ValidError))
-
- # -------- Caclulation of BME as accuracy metric -------
- # Check if data is provided
- if len(ObservationData) != 0:
- # Calculate the initial BME
- BME, KLD, Posterior, DistHellinger = self.BME_Calculator(metamodel, ObservationData, TotalSigma2)
- print('\n')
- print("Updated BME:", BME)
- print("Updated KLD:", KLD)
- print('\n')
-
- # Plot some snapshots of the posterior
- stepSnapshot = self.ExpDesign.stepSnapshot
- if PostSnapshot and postcnt % stepSnapshot == 0:
- MAP = self.ExpDesign.MAP
- parNames = self.ExpDesign.parNames
- print('Posterior snapshot is being plotted...')
- self.posteriorPlot(Posterior, MAP, parNames, 'SeqPosterior_'+str(postcnt))
- postcnt += 1
-
- # Store the updated BME & KLD &
- # Check if data is provided
- if len(ObservationData) != 0:
- SeqBME = np.vstack((SeqBME, BME))
- SeqKLD = np.vstack((SeqKLD, KLD))
- SeqDistHellinger= np.vstack((SeqDistHellinger, DistHellinger))
-
- itrNr += 1
- # Store updated ModifiedLOO and BME in dictonary
- strKey = UtilityFunction + '_rep_' +str(repIdx+1)
- if len(self.validModelRuns) != 0:
- metamodel.seqValidError[strKey] = SeqValidError
-
- # Check if data is provided
- if len(ObservationData) != 0:
- metamodel.SeqBME[strKey] = SeqBME
- metamodel.SeqKLD[strKey] = SeqKLD
- if len(self.validLikelihoods) != 0:
- metamodel.SeqDistHellinger[strKey] = SeqDistHellinger
-
- return metamodel
- #--------------------------------------------------------------------------------------------------------
-
- def create_metamodel(self, Model):
-
- self.ModelObj = Model
- self.ExpDesign.initNrSamples = self.ExpDesign.NrSamples
-
- if self.ExpDesign.Method == 'sequential':
- if self.metaModel.lower() == 'gpe':
- metamodel = self.create_GPE_SeqDesign(Model)
- else:
- metamodel = self.create_PCE_SeqDesign(Model)
-
- elif self.ExpDesign.Method == 'normal':
- metamodel = self.create_PCE_normDesign(Model)
-
- else:
-
- print("The method for experimental design you requested has not been\
- implemented.")
-
- # Cleanup
- # Zip the subdirectories
- try:
- dir_name = Model.Name
- key = Model.Name + '_'
- Model.zip_subdirs(dir_name, key)
- except:
- pass
-
-
- if not self.slicingforValidation:
- return metamodel
-
- # Slice if slicingforValidation is set True
- elif self.slicingforValidation:
-
- # ---- Slice ------
- # Slice the dict based on the index given
- self.index = self.index if type(self.index) is list else [self.index]*len(Model.Output.Names)
-
- OutputDictCalib = {key: self.ModelOutputDict[key][:self.index[idx]] for idx, key in enumerate([list(self.ModelOutputDict.keys())[0]])}
- OutputDictValid = {key: self.ModelOutputDict[key][self.index[idx]:] for idx, key in enumerate([list(self.ModelOutputDict.keys())[0]])}
-
- DictCalib = {key: self.ModelOutputDict[key][:,:self.index[idx]] for idx, key in enumerate(Model.Output.Names)}
- DictValid = {key: self.ModelOutputDict[key][:,self.index[idx]:] for idx, key in enumerate(Model.Output.Names)}
-
- OutputDictCalib.update(DictCalib)
- OutputDictValid.update(DictValid)
-
- # ---- Calibration ----
- PCEModelCalib = self.Slice(metamodel, OutputDictCalib)
-
- # ---- Validation ----
- PCEModelValid = self.Slice(metamodel, OutputDictValid, self.index)
-
-
- return PCEModelCalib , PCEModelValid
-
-
-
- #--------------------------------------------------------------------------------------------------------
- def Slice(self, PCEModel, OutputDict, index_list=None):
- # Deep copy object
- from copy import deepcopy
- newPCEModel = deepcopy(PCEModel)
-
- # Initilize
- newPCEModel.ModelOutputDict = OutputDict
- newPCEModel.CoeffsDict = self.AutoVivification()
- newPCEModel.BasisDict = self.AutoVivification()
- newPCEModel.ScoresDict = self.AutoVivification()
- newPCEModel.clf_poly = self.AutoVivification()
-
- for i, key in enumerate(list(OutputDict.keys())[1:]):
- OutputMatrix = OutputDict[key]
- index = 0 if index_list is None else index_list[i]
- 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)]
- newPCEModel.ScoresDict[key]["y_"+str(index+idx+1)] = PCEModel.ScoresDict[key]["y_"+str(index+idx+1)]
- newPCEModel.clf_poly[key]["y_"+str(index+idx+1)] = PCEModel.clf_poly[key]["y_"+str(index+idx+1)]
- return newPCEModel