diff --git a/BayesValidRox/surrogate_models/ExpDesigns.py b/BayesValidRox/surrogate_models/ExpDesigns.py
index 95b05ebac631dc65cf14579a435b51f6ae69ad8e..4cb853f86e0cdc92b34cb5e238d6bd78da8a8d93 100644
--- a/BayesValidRox/surrogate_models/ExpDesigns.py
+++ b/BayesValidRox/surrogate_models/ExpDesigns.py
@@ -77,15 +77,14 @@ class ExpDesigns:
## Sample the distribution of parameters
if len(Inputs.Marginals[0].InputValues) == 0 and self.SamplingMethod != 'user':
- ## Case I
- # polytype not arbitrary
+ ## Case I = polytype not arbitrary
# Create ExpDesign in the actual space
- X = chaospy.generate_samples(NrSamples, domain=self.JDist, rule=SamplingMethod).T
+ X = chaospy.generate_samples(NrSamples, domain=self.JDist , rule=SamplingMethod).T
else:
- # Case II
+ # Case II:
# polytype arbitrary or Input values are directly prescribed by the user.
# Generate the samples based on requested method
@@ -112,6 +111,7 @@ class ExpDesigns:
origJoints = []
Polytypes = []
+
for parIdx in range(NofPa):
DistType = Inputs.Marginals[parIdx].DistType
params = Inputs.Marginals[parIdx].Parameters
@@ -119,35 +119,40 @@ class ExpDesigns:
if DistType == 'unif':
polytype = 'legendre'
Dist = chaospy.Uniform(lower=params[0],upper=params[1])
-
+
elif DistType == 'norm':
polytype = 'hermite'
Dist = chaospy.Normal(mu=params[0],sigma=params[1])
-
+
elif DistType == 'gamma':
polytype = 'laguerre'
Dist = chaospy.Gamma(shape=params[0],scale=params[1],shift=params[2])
-
+
elif DistType == 'beta':
polytype = 'jacobi'
Dist = chaospy.Beta(alpha=params[0],beta=params[1],
lower=params[2], upper=params[3])
-
+
elif DistType == 'lognormal':
- polytype = 'arbitrary'
- Dist = chaospy.LogNormal(mu=params[0],sigma=params[1])
+ polytype = 'hermite'
+
+ Mu = np.log(params[0]**2 / np.sqrt(params[0]**2 + params[1]**2))
+ Sigma = np.sqrt(np.log(1 + params[1]**2 / params[0]**2))
+
+ Dist = chaospy.LogNormal(mu=Mu,sigma=Sigma)
elif DistType == 'exponential':
polytype = 'arbitrary'
Dist = chaospy.Exponential(scale=params[0],shift=params[1])
-
+
elif DistType == 'weibull':
polytype = 'arbitrary'
Dist = chaospy.Weibull(shape=params[0],scale=params[1],shift=params[2])
-
+
elif DistType is None:
polytype = 'arbitrary'
Dist = None
+
else:
raise ValueError('DistType %s for parameter %s is not available.'%(DistType,parIdx+1))
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 94f6303cc1a1fb03c3709e7bdcec6b9692f02f18..c89d06f327c9506c718ee393844db1cc28b4ab2f 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -178,7 +178,7 @@ class aPCE:
if parIdx_or_parms is None:
print('%s polynomials are parametrically defined! \
Please provide them as an input argument.',polytype)
- parms = parIdx_or_parms
+ Parms = parIdx_or_parms
## Arbitrary polynomials
@@ -202,21 +202,30 @@ class aPCE:
sqrt_b0 = 1
sqrt_bn = lambda n : np.sqrt(n)
bounds = [-np.inf,np.inf]
-
+ polyType = 'Normal'
+ params = [0, 1]
+ Dist = chaospy.Normal(mu=params[0], sigma=params[1])
+
# Uniform distribution
elif polytype == 'legendre':
an = lambda n : np.zeros((n+1))
sqrt_b0 = 1
sqrt_bn = lambda n : np.sqrt(1./(4-(1/np.square(n))))
bounds = [-1,1]
-
+ polyType = 'Uniform'
+ params = [-1, 1]
+ Dist = chaospy.Uniform(params[0], params[1])
+
# Gamma distribution
elif polytype == 'laguerre':
- an = lambda n : 2*n + parms[1]
+ an = lambda n : 2*n + Parms[1]
sqrt_b0 = 1
- sqrt_bn = lambda n : -1*np.sqrt(np.multiply(n ,(n+parms[1]-1)))
+ sqrt_bn = lambda n : -1*np.sqrt(np.multiply(n ,(n+Parms[1]-1)))
bounds = [0,np.inf]
-
+ polyType = 'Gamma'
+ params = self.Inputs.Marginals[parIdx_or_parms].Parameters
+ Dist = chaospy.Gamma(shape=params[0],scale=params[1],shift=params[2])
+
# elif polytype == 'jacobi': # Beta distribution
# # in order to avoid zeros on the denominator some special
# # cancelation cases are hard-coded.
@@ -243,6 +252,8 @@ class aPCE:
# sqrt_bn = @(n) sqrt( 4 .* n .* (n+a) .* (n+b) .* (n+bpa)./((2 .* n + bpa ).^2 .* (2.*n + bpa + 1).*(2.*n + bpa -1)) ) .* 0.5;
#
# bounds = [0,1];
+ # polyType = 'Beta'
+ # prams = [-1, 1]
elif polytype == 'zero': #case 'zero'
# the very special case of constant - recurrence terms are zero:
@@ -250,15 +261,16 @@ class aPCE:
an = lambda n: np.zeros((1, len(n))) #@(n) zeros(1,length(n));
sqrt_bn = lambda n: np.zeros((1, len(n))) #@(n) zeros(1,length(n));
bounds = [-np.inf,np.inf]
-
+ polyType = 'Constant'
+ prams = [-1, 1]
+
else:
raise ValueError('Unknown polynomial type!')
## Assemble the recurrence coefficients into the output
- AB['recurrence terms'] = np.vstack((an(n_max) , np.hstack((sqrt_b0, sqrt_bn(range(1,n_max+1)))))).T
- AB['bounds'] = bounds
+ AB = np.vstack((an(n_max) , np.hstack((sqrt_b0, sqrt_bn(range(1,n_max+1)))))).T
- return AB
+ return AB, Dist
def eval_rec_rule(self,normX,AB,nonrecursive=False):
"""
@@ -337,20 +349,18 @@ class aPCE:
# ----------------
else: # Only for Uniform, Normal, Gamma, Beta
-
- # Forward Rosenblatt transformation
- u = origSpaceDist.fwd(ExpDesignX)
for parIdx in range(NofPa):
- AB = self.poly_rec_coeffs(n_max, polytypes[parIdx], parIdx)
- bounds = AB['bounds']
+ AB, dist = self.poly_rec_coeffs(n_max, polytypes[parIdx], parIdx)
+
+ # Forward Rosenblatt transformation
+ u = origSpaceDist[parIdx].fwd(ExpDesignX[:,parIdx])
# Rosenblatt inverse Transformation based the bounds given by poly_rec_coeffs
- dist = chaospy.Uniform(bounds[0], bounds[1])
- normX = dist.inv(u[:,parIdx])
+ U = dist.inv(u)
- univ_vals[:,parIdx,:] = self.eval_rec_rule(normX,AB['recurrence terms'],nonrecursive=False)
-
+ univ_vals[:,parIdx,:] = self.eval_rec_rule(U,AB,nonrecursive=False)
+
return univ_vals
#--------------------------------------------------------------------------------------------------------
def PCE_create_Psi(self,BasisIndices,univ_p_val):
@@ -777,12 +787,11 @@ class aPCE:
UtilMethod = var
OutputDictY = self.ExpDesign.Y
OutputNames = list(OutputDictY.keys())[1:]
- print("UtilMethod:", UtilMethod)
+
if UtilMethod == 'Entropy':
# ----- Entropy/MMSE -----
# Compute perdiction variance of the old model
Y_PC_can, std_PC_can = self.eval_PCEmodel(X_can)
- print("Y_PC_can:\n", Y_PC_can)
canPredVar = {key:std_PC_can[key]**2 for key in OutputNames}
@@ -941,7 +950,6 @@ class aPCE:
def Utility_runner(self, method, Model, allCandidates,index, sigma2Dict=None, var=None):
if method == 'VarOptDesign':
- print("\nNumber of samples to explore:", len(allCandidates))
U_J_d = self.util_VarBasedDesign(allCandidates, index, var)
elif method == 'BayesOptDesign':
@@ -1229,6 +1237,7 @@ class aPCE:
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]
@@ -1429,10 +1438,8 @@ class aPCE:
# Perdiction
try: # with error bar
clf_poly = self.clf_poly[Outkey][Inkey]
- print("clf_poly:", clf_poly)
- print("PSI_Val:\n",PSI_Val)
y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
- print("y_mean:\n",y_mean)
+
PCEOutputs_mean[:, idx] = y_mean
PCEOutputs_std[:, idx] = y_std
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
index 7949a8f3d96747eae1c17dcd43db06cdf7bb253f..7b0386265f08c3d7d10b467c359ca53460cc175d 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
@@ -100,7 +100,7 @@ if __name__ == "__main__":
MetaModelOpts.RegMethod = 'BRR'
# Print summary of the regression results
- MetaModelOpts.DisplayFlag = True
+ #MetaModelOpts.DisplayFlag = True
# ------ Experimental Design --------
# Generate an experimental design of size NrExpDesign based on a latin
@@ -118,7 +118,7 @@ if __name__ == "__main__":
MetaModelOpts.ExpDesign.SamplingMethod = 'halton'
# Sequential experimental design (needed only for sequential ExpDesign)
- MetaModelOpts.ExpDesign.NrofNewSample = 1
+ MetaModelOpts.ExpDesign.NrofNewSample = 2
MetaModelOpts.ExpDesign.MaxNSamples = 50
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-6
diff --git a/BayesValidRox/tests/BeamTest/SSBeam_Deflection.inp b/BayesValidRox/tests/BeamTest/SSBeam_Deflection.inp
index d9949aca2bee48780b84aabcfecb66f3cce72ecf..00b170341d3c547ee6c7e4f2ed881ab6c2b215da 100644
--- a/BayesValidRox/tests/BeamTest/SSBeam_Deflection.inp
+++ b/BayesValidRox/tests/BeamTest/SSBeam_Deflection.inp
@@ -2,5 +2,5 @@
0.15 % b in m
0.3 % h in m
5 % L in m
-30 % E in Pa
-10 % p in N/m
+30000e+6 % E in Pa
+10000 % p in N/m
diff --git a/BayesValidRox/tests/BeamTest/Test_Beam.py b/BayesValidRox/tests/BeamTest/Test_Beam.py
index c78313df3adc04e16495313b3b9bf947f2a2155b..c3e254fbe763f0de4b568b2665b2d1ac174ef502 100644
--- a/BayesValidRox/tests/BeamTest/Test_Beam.py
+++ b/BayesValidRox/tests/BeamTest/Test_Beam.py
@@ -62,23 +62,23 @@ if __name__ == "__main__":
Inputs.addMarginals()
Inputs.Marginals[0].Name = 'Beam width'
- Inputs.Marginals[0].DistType = 'lognorm'
- Inputs.Marginals[0].Moments = [0.15, 0.0075]
+ Inputs.Marginals[0].DistType = 'lognormal'
+ Inputs.Marginals[0].Parameters = [0.15, 0.0075]
Inputs.addMarginals()
Inputs.Marginals[1].Name = 'Beam height'
- Inputs.Marginals[1].DistType = 'lognorm'
- Inputs.Marginals[1].Moments = [0.3, 0.015]
+ Inputs.Marginals[1].DistType = 'lognormal'
+ Inputs.Marginals[1].Parameters = [0.3, 0.015]
Inputs.addMarginals()
Inputs.Marginals[2].Name = 'Youngs modulus'
- Inputs.Marginals[2].DistType = 'norm' #'lognorm'
- Inputs.Marginals[2].Moments = [30, 4.5] #[30000e+6, 4500e+6]
+ Inputs.Marginals[2].DistType = 'lognormal'
+ Inputs.Marginals[2].Parameters = [30000e+6, 4500e+6]
Inputs.addMarginals()
Inputs.Marginals[3].Name = 'Uniform load'
- Inputs.Marginals[3].DistType = 'norm' #'lognorm'
- Inputs.Marginals[3].Moments = [10, 2]
+ Inputs.Marginals[3].DistType = 'lognormal'
+ Inputs.Marginals[3].Parameters = [1e4, 2e3]
#=====================================================
#====== POLYNOMIAL CHAOS EXPANSION METAMODELS ======
@@ -107,49 +107,68 @@ 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'
MetaModelOpts.ExpDesign.NrSamples = 20
- MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
- MetaModelOpts.ExpDesign.Method = 'sequential' # 1) normal 2) sequential
- #MetaModelOpts.ExpDesign.X = np.load('CollocationPoints.npy')
+
+ # 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 = 'halton'
+
+ #MetaModelOpts.ExpDesign.X = np.load('EDX_Beam9points.npy')
# Sequential experimental design (needed only for sequential ExpDesign)
- MetaModelOpts.ExpDesign.MaxNSamples = 50 #150
- MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-3
+ MetaModelOpts.ExpDesign.NrofNewSample = 2
+ MetaModelOpts.ExpDesign.MaxNSamples = 50
+ MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-6
+
+ DiscrepancyOpts = Discrepancy('')
+ DiscrepancyOpts.Type = 'Gaussian'
+ DiscrepancyOpts.Parameters = np.square([1e-8, 1e-2, 1e-2, 1e-3, 4e-3, 5e-3, 1e-3,
+ 1e-2, 1e-2, 1e-2, 1e-2, 1e-9])
+ MetaModelOpts.Discrepancy = DiscrepancyOpts
# Plot the posterior snapshots for SeqDesign
MetaModelOpts.ExpDesign.PostSnapshot = True
MetaModelOpts.ExpDesign.stepSnapshot = 1
- MetaModelOpts.ExpDesign.MAP = (0.150064, 0.299698, 30.763206, 10.164872)
+ MetaModelOpts.ExpDesign.MAP = (0.150064, 0.299698, 30.763206e9, 10.164872e2)
MetaModelOpts.ExpDesign.parNames = ['Beam width', 'Beam height', 'Youngs modulus', 'Uniform load']
-
- # -------- Optimality criteria: Optimization ------
- # 1)'dual annealing' 2)'minimization' 3)'BayesOptDesign'
- MetaModelOpts.ExpDesign.SeqOptimMethod = 'BayesOptDesign'
- MetaModelOpts.ExpDesign.ExplorationMethod = 'Voronoi' #1)'Voronoi' 2)'LHS' 3) 'MC'
+
+ # ------------------------------------------------
+ # ------- Sequential Design configuration --------
+ # ------------------------------------------------
+ # 1) 'None' 2) 'epsilon-decreasing'
+ MetaModelOpts.ExpDesign.TradeOffScheme = 'None'
+ #MetaModelOpts.ExpDesign.nReprications = 2
+ # -------- Exploration ------
+ #1)'Voronoi' 2)'MC' 3)'LHS' 4)'dual annealing'
+ MetaModelOpts.ExpDesign.ExploreMethod = 'Voronoi'
+
+ # Use when 'dual annealing' chosen
+ MetaModelOpts.ExpDesign.MaxFunItr = 200
+
+ # Use when 'Voronoi' or 'MC' or 'LHS' chosen
+ MetaModelOpts.ExpDesign.NCandidate = 1000
+ MetaModelOpts.ExpDesign.NrofCandGroups = 4
- MetaModelOpts.ExpDesign.MaxFunItr = 100
+ # -------- Exploitation ------
+ # 1)'BayesOptDesign' 2)'VarOptDesign' 3)'alphabetic' 4)'Space-filling'
+ MetaModelOpts.ExpDesign.ExploitMethod = 'VarOptDesign'
- MetaModelOpts.ExpDesign.NCandidate = 1000 # 5000
- MetaModelOpts.ExpDesign.NrofCandGroups = 4
+ # BayesOptDesign -> when data is available
# 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
# 3)APP (A-Posterior-percision)
- MetaModelOpts.ExpDesign.UtilityFunction = 'DKL'
-
-# DiscrepancyOpts = Discrepancy('')
-# DiscrepancyOpts.Type = 'Gaussian'
-# DiscrepancyOpts.Parameters = [1e-8, 1e-2, 1e-2, 1e-3, 4e-3, 5e-3, 1e-3,
-# 1e-2, 1e-2, 1e-2, 1e-2, 1e-9] #1e-5
-# MetaModelOpts.Discrepancy = DiscrepancyOpts
+ #MetaModelOpts.ExpDesign.UtilityFunction = 'DKL' #['DKL', 'DPP']
- # -------- Optimality criteria: alphabetic ------
-# MetaModelOpts.ExpDesign.SeqOptimMethod = 'alphabetic'
-# MetaModelOpts.ExpDesign.NCandidate = 5000
-#
-# # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
-# # 3)K-Opt (K-Optimality)
-# MetaModelOpts.ExpDesign.UtilityFunction = 'D-Opt'
-#
+ # VarBasedOptDesign -> when data is not available
+ # Only with Vornoi >>> 1)Entropy 2)EIGF, 3)ALM, 4)LOOCV
+ MetaModelOpts.ExpDesign.UtilityFunction = 'Entropy'#['EIGF', 'Entropy', 'LOOCV']
+ # alphabetic
+ # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
+ # 3)K-Opt (K-Optimality)
+ #MetaModelOpts.ExpDesign.UtilityFunction = 'D-Opt' #['D-Opt', 'A-Opt', 'K-Opt']
# >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Adaptive sparse arbitrary polynomial chaos expansion
diff --git a/BayesValidRox/tests/BeamTest/myBeam9points b/BayesValidRox/tests/BeamTest/myBeam9points
index ff025184b64135d4af02fdbe06ff57e9f6251e4b..559357bf494373c8ebc0a02892f79831dd70313d 100755
Binary files a/BayesValidRox/tests/BeamTest/myBeam9points and b/BayesValidRox/tests/BeamTest/myBeam9points differ
diff --git a/BayesValidRox/tests/BeamTest/myBeam9points.cpp b/BayesValidRox/tests/BeamTest/myBeam9points.cpp
index 2ad6b9e6ee65dec6016fe3bafecf8f146932d7f6..546428bd69ea0243cab21ab8f4ca3e59de4d7a3c 100644
--- a/BayesValidRox/tests/BeamTest/myBeam9points.cpp
+++ b/BayesValidRox/tests/BeamTest/myBeam9points.cpp
@@ -66,8 +66,8 @@ int main(int argc, char* argv[]) {
beam_width = params[0];
beam_height = params[1];
beam_span = params[2];
- youngs_modulus = params[3] * 1e+9;
- load = params[4] * 1e+3;
+ youngs_modulus = params[3];
+ load = params[4];
inputfile.close();
} else {