diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 76aede9cc06f5d8228c4658fe4d018d72c3e9c95..b95dd59e2848835c2ccdc6ac3fb3c74dbe70c556 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -21,6 +21,8 @@ plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=24)
plt.rc('ytick', labelsize=24)
plt.rc('axes', labelsize=24)
+plt.rc('axes', grid=True)
+plt.rc('grid', linestyle=":")
plt.rc('savefig', dpi=1000)
diff --git a/BayesValidRox/PyLink/PyLinkForwardModel.py b/BayesValidRox/PyLink/PyLinkForwardModel.py
index b91e37fef8852bd4ecf6f4fda64b13ef7117bad7..9a3f166f7d65b002320751abdb94dcd5d91ea988 100644
--- a/BayesValidRox/PyLink/PyLinkForwardModel.py
+++ b/BayesValidRox/PyLink/PyLinkForwardModel.py
@@ -202,13 +202,9 @@ class PyLinkForwardModel:
#os.rename(self.AuxFile, NewShellfile)
- #--------- Simulation Starts here --------------------
# Update the parametrs in Input file
self.UpdateParamterNew(NewInputfile, CollocationPoints)
- #-----------------------------------------------------
- # Simulation Starts here
- #-----------------------------------------------------
# Update the user defined command and the execution path
try:
NewCommand = self.Command.replace(self.InputFile[0], NewInputfile[0]).replace(self.InputFile[1], NewInputfile[1])
diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index cad82103fd8c1b7a76d28d4d79113438c51d2377..2fdbc405c34bdc632153d4b8069d84c4ea8a516d 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -95,7 +95,10 @@ class RegressionFastARD():
----------
n_iter: int, optional (DEFAULT = 100)
Maximum number of iterations
-
+
+ start: list, optional (DEFAULT = None)
+ Initial selected features.
+
tol: float, optional (DEFAULT = 1e-3)
If absolute change in precision parameter for weights is below threshold
algorithm terminates.
@@ -147,9 +150,10 @@ class RegressionFastARD():
'''
- def __init__( self, n_iter = 300, tol = 1e-3, fit_intercept = True,
+ def __init__( self, n_iter = 300, start=None, tol = 1e-3, fit_intercept = True,
copy_X = True, compute_score=False, verbose = False):
self.n_iter = n_iter
+ self.start = start
self.tol = tol
self.scores_ = list()
self.fit_intercept = fit_intercept
@@ -213,18 +217,28 @@ class RegressionFastARD():
A = np.PINF * np.ones(n_features)
active = np.zeros(n_features , dtype = np.bool)
- # in case of almost perfect multicollinearity between some features
- # start from feature 0
- if np.sum( XXd - X_mean**2 < np.finfo(np.float32).eps ) >= 0:
- A[0] = np.finfo(np.float16).eps
- active[0] = True
- else:
- # start from a single basis vector with largest projection on targets
+
+ if self.start is not None:
+ start = self.start
+ # start from a given start basis vector
proj = XY**2 / XXd
- start = np.argmax(proj)
active[start] = True
A[start] = XXd[start]/( proj[start] - var_y)
-
+
+ else:
+ # in case of almost perfect multicollinearity between some features
+ # start from feature 0
+ if np.sum( XXd - X_mean**2 < np.finfo(np.float32).eps ) >= 0:
+ A[0] = np.finfo(np.float16).eps
+ active[0] = True
+
+ else:
+ # start from a single basis vector with largest projection on targets
+ proj = XY**2 / XXd
+ start = np.argmax(proj)
+ active[start] = True
+ A[start] = XXd[start]/( proj[start] - var_y)
+
warning_flag = 0
scores_ = []
for i in range(self.n_iter):
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index b4eeb1efd2d34486ed4e39fc16d964fed927cb06..18ee23cb2fb5a94c40c7a841af5ff5e78f968ae3 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -25,6 +25,8 @@ plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=24)
plt.rc('ytick', labelsize=24)
plt.rc('axes', labelsize=24)
+plt.rc('axes', grid=True)
+plt.rc('grid', linestyle="-")
plt.rc('savefig', dpi=1000)
import scipy as sp
@@ -34,6 +36,7 @@ from itertools import combinations, product
from scipy import stats
from sklearn.linear_model import LinearRegression, ARDRegression
from sklearn import linear_model
+from sklearn.metrics import r2_score
import os
import sys
import seaborn as sns
@@ -367,8 +370,8 @@ class aPCE:
origSpaceDist = self.ExpDesign.JDist
if ExpDesignX.ndim !=2: ExpDesignX = ExpDesignX.reshape(1, len(ExpDesignX))
NofSamples, NofPa = ExpDesignX.shape
- P = self.MaxPceDegree + 1
n_max = self.MaxPceDegree
+ P = n_max + 1
# Allocate the output matrix
univ_vals = np.zeros((NofSamples,NofPa, P))
@@ -435,7 +438,7 @@ class aPCE:
Psi[:,aa] = np.multiply(Psi[:,aa] , np.reshape(univ_p_val[:, m, basisIdx], Psi[:,aa].shape))
except:
raise ValueError
-
+
return Psi
#--------------------------------------------------------------------------------------------------------
@@ -535,9 +538,9 @@ class aPCE:
- # Adaptive aPCE (not Bayes sparse)adaptive aPCE
+ # Ordinary least square method (OLS)
else:
- PolynomialDegrees_Sparse = BasisIndices
+ PolynomialDegrees_Sparse = BasisIndices.toarray() if sparsity else BasisIndices
PSI_Sparse = PSI
#==============================================================================
#==== Computation of the expansion coefficients for arbitrary Polynomial Chaos
@@ -559,14 +562,14 @@ class aPCE:
# Initialization
- self.DegreeArray = np.array(range(self.MinPceDegree,self.MaxPceDegree+1))
- qnorm = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
qAllCoeffs, AllCoeffs = {}, {}
qAllIndices_Sparse ,AllIndices_Sparse = {}, {}
qAllclf_poly, Allclf_poly = {}, {}
qAllnTerms, AllnTerms = {}, {}
qAllLCerror, AllLCerror = {}, {}
+ qnorm = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
+
errorIncreases = False
# Overall score indicator: this will be used to compare between different
# methods
@@ -597,8 +600,9 @@ class aPCE:
for degIdx, deg in enumerate(self.DegreeArray):
for qidx, q in enumerate(qnorm):
- # Generate the polynomial basis indices
- BasisIndices = self.PolyBasisIndices(degree = deg, q=q)
+
+ # 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,univ_p_val)
@@ -707,7 +711,6 @@ class aPCE:
#--------------------------------------------------------------------------------------------------------
def CoeffsUpdateaPCE(self, CollocationPoints, ModelOutput, BasisIndices):
-
# Evaluate the univariate polynomials on ExpDesig
univ_p_val = self.univ_basis_vals(CollocationPoints)
@@ -715,7 +718,7 @@ class aPCE:
Psi = self.PCE_create_Psi(BasisIndices,univ_p_val)
- PolynomialDegrees, PSI_Sparse, Coeffs, clf_poly = self.RS_Builder(Psi, ModelOutput, BasisIndices)
+ PolynomialDegrees, PSI_Sparse, Coeffs, clf_poly = self.RS_Builder(Psi, ModelOutput, BasisIndices) #,startIdxes=startIdxes)
# Calculate and save the score of LOOCV
LOOCVScore, LCerror = self.CorrectedLeaveoutError(PSI_Sparse, Coeffs, ModelOutput, clf_poly)
@@ -871,11 +874,21 @@ class aPCE:
self.ScoresDict = self.AutoVivification()
self.clf_poly = self.AutoVivification()
self.LCerror = self.AutoVivification()
+ self.allBasisIndices = self.AutoVivification()
# Read observations or MCReference
if Model.MeasurementFile is not None or len(Model.Observations) != 0:
self.Observations = Model.read_Observation()
-
+
+ # Generate Basis indices matrix
+ self.DegreeArray = np.array(range(self.MinPceDegree,self.MaxPceDegree+1))
+ qnorm = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
+
+ for degIdx, deg in enumerate(self.DegreeArray):
+ for qidx, q in enumerate(qnorm):
+ # Generate the polynomial basis indices
+ self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree = deg, q=q)
+
# ---- Experimental Design ------
# Get the collocation points to run the forward model
@@ -2083,7 +2096,7 @@ class aPCE:
#infEntropy = logBME - postExpPrior - postExpLikelihoods
# If PostSnapshot is True, plot likelihood vs refrence
- if PostSnapshot or len(validLikelihoods)!=0:
+ if PostSnapshot and 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')
@@ -2119,16 +2132,13 @@ class aPCE:
# Run the PCE model with the generated samples
PCEOutputs, PCEOutputs_std = self.eval_PCEmodel(Samples)
- validError_dict = self.AutoVivification()
+ validError_dict = {}
# Loop over the keys and compute RMSE error.
for key in OutputName:
- # Residuals
- residual = PCEOutputs[key] - ModelOutputs
- validErr = np.mean(residual**2, axis=0)/np.var(ModelOutputs, axis=0, ddof=1)
+ validErr = 1 - r2_score(PCEOutputs[key],ModelOutputs,multioutput='variance_weighted')
- for idx in range(ModelOutputs.shape[1]):
- validError_dict[key]["y_"+str(idx+1)] = validErr[idx]
-
+ validError_dict[key] = validErr
+
return validError_dict
#--------------------------------------------------------------------------------------------------------
@@ -2168,6 +2178,7 @@ class aPCE:
return RMSE_Mean, RMSE_std
#--------------------------------------------------------------------------------------------------------
def create_PCE_SeqDesign(self, Model):
+
#Initialization
errorIncreases = False
@@ -2184,7 +2195,7 @@ class aPCE:
if not isinstance(self.ExpDesign.UtilityFunction, list):
UtilityFunctions = [self.ExpDesign.UtilityFunction]
- # Initial PCEModel
+ # ---------- Initial PCEModel ----------
PCEModel = self.create_PCE_normDesign(Model)
initPCEModel= PCEModel
@@ -2227,19 +2238,19 @@ class aPCE:
initYprev = deepcopy(initPCEModel.ModelOutputDict)
# Read the initial ModifiedLOO
- Scores_all=[]
+ Scores_all, varExpDesignY =[], []
for OutName in OutputName:
Scores_all.append(list(initPCEModel.ScoresDict[OutName].values()))
+ varExpDesignY.append(np.var(initPCEModel.ExpDesign.Y[OutName],axis=0))
Scores = [item for sublist in Scores_all for item in sublist]
- initModifiedLOO = [1-score for score in Scores]
+ 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_()
- Scores_all=[]
- for OutName in OutputName:
- Scores_all.append(list(initValidError[OutName].values()))
- initValidError = [item for sublist in Scores_all for item in sublist]
-
+ initValidError = list(initValidError.values())
+ print("\nInitial ValidError:", initValidError)
+
self.SeqModifiedLOO = {}
self.SeqBME = {}
self.seqRMSEMean = {}
@@ -2330,25 +2341,24 @@ class aPCE:
# -------- SparseBayes calculation of the coefficients-------
prevPCEModel = PCEModel
- PCEModel = PCEModel.create_PCE_normDesign(Model)
-
+ PCEModel = PCEModel.create_PCE_normDesign(Model)#, OptDesignFlag=True)
+
# Compute the validation error
if len(self.validModelRuns) != 0:
validError = PCEModel.validError_()
- Scores_all=[]
- for OutName in OutputName:
- Scores_all.append(list(validError[OutName].values()))
- ValidError = [item for sublist in Scores_all for item in sublist]
+ ValidError = list(validError.values())
print("\nUpdated ValidError:", ValidError)
- print("\n")
# Extract Modified LOO from Output
- Scores_all = []
+ Scores_all, varExpDesignY =[], []
for OutName in OutputName:
Scores_all.append(list(PCEModel.ScoresDict[OutName].values()))
+ varExpDesignY.append(np.var(PCEModel.ExpDesign.Y[OutName],axis=0))
Scores = [item for sublist in Scores_all for item in sublist]
- ModifiedLOO = [1-score for score in Scores]
-
+ 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)
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalyticFunc_Demo.py b/BayesValidRox/tests/AnalyticalFunction/AnalyticFunc_Demo.py
index 4290dd4c83dae1babdc69b7b4b90390b4ef0f26e..f369086923cbe9183c0e0536167a2131c9920bb5 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalyticFunc_Demo.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalyticFunc_Demo.py
@@ -68,9 +68,10 @@ if __name__ == "__main__":
# Select the sparse least-square minimization method for
# the PCE coefficients calculation:
- # 1)AaPCE: Adaptive aPCE 2)BRR: Bayesian Ridge Regression
- # 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression (Not Working)
- MetaModelOpts.RegMethod = 'AaPCE'
+ # 1)OLS: Ordinary Least Square 2)BRR: Bayesian Ridge Regression
+ # 3)LARS: Least angle regression 4)ARD: Bayesian ARD Regression
+ # 5)SGDR: Stochastic gradient descent learning
+ MetaModelOpts.RegMethod = 'OLS'
# Print summary of the regression results
#MetaModelOpts.DisplayFlag = True
@@ -80,10 +81,15 @@ if __name__ == "__main__":
# hypercube sampling of the input model or user-defined values of X and/or Y:
MetaModelOpts.addExpDesign()
+ # One-shot (normal) or Sequential Adaptive (sequential) Design
+ MetaModelOpts.ExpDesign.Method = 'normal'
NrofInitSamples = 300
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
- MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
- MetaModelOpts.ExpDesign.Method = 'normal' # 1) normal 2) sequential
+
+ # Sampling methods
+ # 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) chebyshev(FT)
+ # 7) korobov 8) grid(FT) 9) nested_grid(FT) 10)user
+ MetaModelOpts.ExpDesign.SamplingMethod = 'latin_hypercube'
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.MaxNSamples = 100
@@ -105,26 +111,26 @@ if __name__ == "__main__":
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Adaptive sparse arbitrary polynomial chaos expansion
- PCEModel_AaPCE = MetaModelOpts.create_PCE(Model)
+ PCEModel_OLS = MetaModelOpts.create_PCE(Model)
# Extract the basis and coefficients
- Basis_AaPCE = PCEModel_AaPCE.BasisDict['Z']['y_2']
- COFFS_AaPCE = PCEModel_AaPCE.CoeffsDict['Z']['y_2']
+ Basis_OLS = PCEModel_OLS.BasisDict['Z']['y_2']
+ COFFS_OLS = PCEModel_OLS.CoeffsDict['Z']['y_2']
#=====================================================
#=========== DEMO OF BAYESPCE METAMODEL ===========
#=====================================================
fig1=plt.figure()
barWidth = 0.3
- plt.title("Bar plot of the weights (AaPCE)")
+ plt.title("Bar plot of the weights (OLS)")
# Create cyan bars
- r = np.arange(len(Basis_AaPCE))
- plt.bar(r, COFFS_AaPCE, width = barWidth, color = 'cyan', bottom=0.001,
+ r = np.arange(len(Basis_OLS))
+ plt.bar(r, COFFS_OLS, width = barWidth, color = 'cyan', bottom=0.001,
edgecolor = 'black', capsize=7, label='Cofficients')
# general layout
- plt.xticks([r + barWidth for r in range(len(r))], ['term'+str(i+1) for i in range(len(Basis_AaPCE))])
+ plt.xticks([r + barWidth for r in range(len(r))], ['term'+str(i+1) for i in range(len(Basis_OLS))])
plt.xlabel("Features")
plt.ylabel("Values of the weights")
@@ -152,10 +158,10 @@ if __name__ == "__main__":
# ------ Experimental Design --------
BRRMetaModelOpts.addExpDesign()
- BRRMetaModelOpts.ExpDesign.NrSamples = PCEModel_AaPCE.ExpDesign.X.shape[0]
+ BRRMetaModelOpts.ExpDesign.NrSamples = PCEModel_OLS.ExpDesign.X.shape[0]
BRRMetaModelOpts.ExpDesign.SamplingMethod = 'user'
- BRRMetaModelOpts.ExpDesign.X = PCEModel_AaPCE.ExpDesign.X
- BRRMetaModelOpts.ExpDesign.Y = PCEModel_AaPCE.ExpDesign.Y
+ BRRMetaModelOpts.ExpDesign.X = PCEModel_OLS.ExpDesign.X
+ BRRMetaModelOpts.ExpDesign.Y = PCEModel_OLS.ExpDesign.Y
BRRMetaModelOpts.ExpDesign.Method = 'normal'
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index 84eaa5dd7607e2c25fb5f8b93ed7f86dcca83703..ad37276fa1540534f0f35b4b8c6de8ccfb103544 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -9,14 +9,14 @@ Created on Fri Aug 9 16:25:43 2019
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
-plt.style.use('seaborn-white')
-plt.rcParams.update({'font.size': 18})
+# plt.style.use('seaborn-paper')
+plt.rcParams.update({'font.size': 24})
plt.rc('font', family='sans-serif', serif='Arial')
plt.rc('figure', figsize = (24, 16))
plt.rc('text', usetex=True)
-plt.rc('xtick', labelsize=18)
-plt.rc('ytick', labelsize=18)
-plt.rc('axes', labelsize=18)
+plt.rc('xtick', labelsize=24)
+plt.rc('ytick', labelsize=24)
+plt.rc('axes', labelsize=24)
import seaborn as sns
import os, sys
@@ -39,7 +39,7 @@ from BayesInference.BayesInference import BayesInference, Discrepancy
if __name__ == "__main__":
# Number of parameters
- ndim = 2
+ ndim = 10 # 2, 10
#=====================================================
#============= COMPUTATIONAL MODEL ================
@@ -110,7 +110,7 @@ if __name__ == "__main__":
# One-shot (normal) or Sequential Adaptive (sequential) Design
MetaModelOpts.ExpDesign.Method = 'sequential'
- NrofInitSamples = 10
+ NrofInitSamples = 50
MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
# Sampling methods
@@ -120,7 +120,7 @@ if __name__ == "__main__":
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
- MetaModelOpts.ExpDesign.MaxNSamples = 50
+ MetaModelOpts.ExpDesign.MaxNSamples = 100
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-16
# Defining the measurement error, if it's known a priori
@@ -144,7 +144,7 @@ if __name__ == "__main__":
# MetaModelOpts.ExpDesign.nReprications = 2
# -------- Exploration ------
#1)'Voronoi' 2)'random' 3)'latin_hypercube' 4)'dual annealing'
- MetaModelOpts.ExpDesign.ExploreMethod = 'Voronoi'
+ MetaModelOpts.ExpDesign.ExploreMethod = 'random'
# Use when 'dual annealing' chosen
MetaModelOpts.ExpDesign.MaxFunItr = 200
@@ -155,16 +155,16 @@ if __name__ == "__main__":
# -------- Exploitation ------
# 1)'BayesOptDesign' 2)'VarOptDesign' 3)'alphabetic' 4)'Space-filling'
- MetaModelOpts.ExpDesign.ExploitMethod = 'VarOptDesign'
+ MetaModelOpts.ExpDesign.ExploitMethod = 'BayesOptDesign'
# BayesOptDesign -> when data is available
# 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
# 3)APP (A-Posterior-percision)
- # MetaModelOpts.ExpDesign.UtilityFunction = 'DKL'# ['DKL', 'BME', 'infEntropy']
+ MetaModelOpts.ExpDesign.UtilityFunction = 'DKL'# ['DKL', 'BME', 'infEntropy']
# VarBasedOptDesign -> when data is not available
# Only with Vornoi >>> 1)Entropy 2)EIGF, 3)ALM, 4)LOOCV
- MetaModelOpts.ExpDesign.UtilityFunction = 'LOOCV'#['EIGF', 'Entropy', 'LOOCV']
+ # MetaModelOpts.ExpDesign.UtilityFunction = 'Entropy'#['EIGF', 'Entropy', 'LOOCV']
# alphabetic
# 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)