diff --git a/BayesValidRox/tests/PA-A/Benchmark_PAA.py b/BayesValidRox/tests/PA-A/Benchmark_PAA.py
index 85028cb719245467c784656d764c9d41e9632659..56c5450b0b4a6a07b6effb069ed2e8f8ef980ca9 100755
--- a/BayesValidRox/tests/PA-A/Benchmark_PAA.py
+++ b/BayesValidRox/tests/PA-A/Benchmark_PAA.py
@@ -94,7 +94,7 @@ def kdePlot_BME(dof, BME, ModelNames, plotName):
plt.savefig('./'+OutputDir+directory+'BME_TOM_'+plotName+'.svg', bbox_inches = 'tight')
plt.savefig('./'+OutputDir+directory+'BME_TOM_'+plotName+'.pdf', bbox_inches = 'tight')
- plt.show()
+
plt.close()
return logTOMBME
@@ -163,10 +163,10 @@ def kdePlot_BayesFactor(BME_Dict, plotName):
ax.fill_between(linearized, 0, Z2, color=Colors[i], alpha=0.25)
#ax.set_ylim([0,0.75])
- # Draw BF thresholds according to Jeffreys 1961
- ax.axvline(x=0,ymin=0, linewidth=4, color='dimgrey') # Strong evidence for both models
- ax.axvline(x=1,ymin=0,linewidth=4, color='orange') # Strong evidence for one model
- ax.axvline(x=2,ymin=0,linewidth=4, color='r') # Decisive evidence for one model
+ # Draw BF significant levels according to Jeffreys 1961
+ ax.axvline(x=np.log10(3),ymin=0, linewidth=4, color='dimgrey') # Strong evidence for both models
+ ax.axvline(x=np.log10(10),ymin=0,linewidth=4, color='orange') # Strong evidence for one model
+ ax.axvline(x=np.log10(100),ymin=0,linewidth=4, color='r') # Decisive evidence for one model
# legend
BF_label = key_i + '/' + key_j
@@ -205,11 +205,11 @@ def kdePlot_BayesFactor(BME_Dict, plotName):
axes[i, 0].set_ylabel('Probability', fontsize=SIZE)
# Adjust subplots
- plt.subplots_adjust(wspace=0.1, hspace=0.1)
+ plt.subplots_adjust(wspace=0.2, hspace=0.1)
plt.savefig('./'+OutputDir+directory+'BayesFactor'+plotName+'.svg', bbox_inches = 'tight')
plt.savefig('./'+OutputDir+directory+'BayesFactor'+plotName+'.pdf', bbox_inches = 'tight')
- plt.show()
+
plt.close()
def BoxPlot_ModelWeights(ModelWeights, plotName):
@@ -226,7 +226,7 @@ def BoxPlot_ModelWeights(ModelWeights, plotName):
filtered_data = [d[m] for d, m in zip(ModelWeights, mask.T)]
# Create the boxplot
- bp = ax.boxplot(filtered_data, patch_artist=True)
+ bp = ax.boxplot(filtered_data, patch_artist=True, showfliers=False)
## change outline color, fill color and linewidth of the boxes
for box in bp['boxes']:
@@ -248,8 +248,8 @@ def BoxPlot_ModelWeights(ModelWeights, plotName):
median.set(color='#b2df8a', linewidth=2)
## change the style of fliers and their fill
- for flier in bp['fliers']:
- flier.set(marker='o', color='#e7298a', alpha=0.75)
+ # for flier in bp['fliers']:
+ # flier.set(marker='o', color='#e7298a', alpha=0.75)
## Custom x-axis labels
ax.set_xticklabels(list(BME_Dict.keys()))
@@ -272,7 +272,6 @@ def BoxPlot_ModelWeights(ModelWeights, plotName):
#fig.savefig('./'+OutputDir+directory+'ModelWeights'+plotName+'.svg', bbox_inches='tight')
fig.savefig('./'+OutputDir+directory+'ModelWeights'+plotName+'.pdf', bbox_inches='tight')
- plt.show()
plt.close()
def perturbData(data,OutputNames,noItr,noiseLevel):
@@ -298,9 +297,9 @@ if __name__ == "__main__":
#==========================================================================
#======================== Set Parameters ==============================
#==========================================================================
- # Set dir name for plots
- OutputDir = (r'Outputs_BayesAnalysis_PAA/')
-
+ # ---------- set scenario params ----------
+ inletLoc = 'left' # left or top
+
# ---------- set BayesValidRox params ----------
PCEExpDesignMethod = 'normal' # 'normal' or 'sequential'
nInitSamples = 500 #50 Initial No. of orig. Model runs for surrogate training
@@ -309,16 +308,16 @@ if __name__ == "__main__":
BootstrapNoise = 0.0025 # Noise amount for bootstraping in Bayesian analysis
# Perturbe data
- data = pd.read_csv('data/stokesDataValid.csv')
+ data = pd.read_csv('data/stokesDataValid_{}Open.csv'.format(inletLoc))
perturbedDataAvg = perturbData(data,['velocity [m/s]', 'p'],nBootstrapItr,BootstrapNoise)
- np.savetxt('./data/perturbedValidDataAvg.csv',perturbedDataAvg, delimiter=',')
- # perturbedDataAvg = np.loadtxt('./data/perturbedValidDataAvg.csv',delimiter=',')
+ np.savetxt('./data/perturbedValidDataAvg_{}Open.csv'.format(inletLoc),perturbedDataAvg, delimiter=',')
+ # perturbedDataAvg = np.loadtxt('./data/perturbedValidDataAvg_{}Open.csv'.format(inletLoc),delimiter=',')
paramsAvg = (nInitSamples,nTotalSamples,nBootstrapItr,perturbedDataAvg)
# Perturbe non-averaged data for PNM
- data = pd.read_csv('data/stokesDataValid_without_averaging.csv')
+ data = pd.read_csv('data/stokesDataValid_{}Open_without_averaging.csv'.format(inletLoc))
perturbedData = perturbData(data,['velocity [m/s]', 'p'],nBootstrapItr,BootstrapNoise)
- np.savetxt('./data/perturbedValidData.csv',perturbedData, delimiter=',')
+ np.savetxt('./data/perturbedValidData_{}Open.csv'.format(inletLoc),perturbedData, delimiter=',')
# perturbedData = np.loadtxt('./data/perturbedValidData.csv',delimiter=',')
params = (nInitSamples,nTotalSamples,nBootstrapItr,perturbedData)
@@ -326,42 +325,48 @@ if __name__ == "__main__":
#==========================================================================
#==================== Run main scripts for PA-B =======================
#==========================================================================
- result_folder = './' #'./Results_22_06_2021'
+ # Set dir name for plots
+ OutputDir = (r'Outputs_BayesAnalysis_PAA_{}/'.format(inletLoc))
+ result_folder = './' #'./Results_20_09_2021'
+
+ # ------------- Run scripts ---------
+ # Stokes-PN model with the averaged data
+ _, _ , BayesValid_PNM = stokespnm.run(paramsAvg, inletLoc,
+ PCEEDMethod=PCEExpDesignMethod)
+
+ # Stokes-PN model without the averaged data
+ _, _, BayesValid_PNM_NA = stokespnm.run(params, averaging=False,inletLoc=inletLoc,
+ PCEEDMethod=PCEExpDesignMethod)
- # Model stokesdarcy-pnm with the averaged data
- # PCEModel_PNM, BayesCalib_PNM, BayesValid_PNM = stokespnm.run(paramsAvg, PCEEDMethod=PCEExpDesignMethod)
- # Load the objects
+ # StokesDarcy with Classical IC (Beaver-Joseph)
+ _, _, BayesValid_BJ = stokesdarcy.run(paramsAvg, couplingcond='BJ',inletLoc=inletLoc,
+ PCEEDMethod=PCEExpDesignMethod)
+
+
+ # StokesDarcy with Generalized (ER) IC
+ _, _, BayesValid_ER = stokesdarcy.run(paramsAvg, couplingcond='ER', inletLoc=inletLoc,
+ PCEEDMethod=PCEExpDesignMethod)
+
+
+ # ------------- Load the objects ---------
+ # Stokes-PN model with the averaged data
# with open(result_folder+'/outputs_ffpm-stokespnm/PA_A_Bayesffpm-stokespnm-valid.pkl', 'rb') as input:
# BayesValid_PNM = joblib.load(input)
- # Model stokesdarcy-pnm without the averaged data
- # PCEModel_PNM_NA, BayesCalib_PNM_NA, BayesValid_PNM_NA = stokespnm.run(params, averaging=False,
- # PCEEDMethod=PCEExpDesignMethod)
-
- # Load the objects
+ # Stokes-PN model without the averaged data
# with open(result_folder+'/outputs_ffpm-stokespnmNA/PA_A_Bayesffpm-stokespnmNA-valid.pkl', 'rb') as input:
# BayesValid_PNM_NA = joblib.load(input)
-
- # Model stokesdarcy-BJ
- PCEModel_BJ, BayesCalib_BJ, BayesValid_BJ = stokesdarcy.run(paramsAvg, couplingcond='BJ',
- PCEEDMethod=PCEExpDesignMethod)
-
- # Load the objects
+
+ # StokesDarcy with Classical IC (Beaver-Joseph)
# with open(result_folder+'/outputs_ffpm-stokesdarcyBJ/PA_A_Bayesffpm-stokesdarcyBJ-valid.pkl', 'rb') as input:
# BayesValid_BJ = joblib.load(input)
-
- # Model stokesdarcy-ER
- PCEModel_ER, BayesCalib_ER, BayesValid_ER = stokesdarcy.run(paramsAvg, couplingcond='ER',
- PCEEDMethod=PCEExpDesignMethod)
-
-
- # Load the objects
+
+ # StokesDarcy with Generalized (ER) IC
# with open(result_folder+'/outputs_ffpm-stokesdarcyER/PA_A_Bayesffpm-stokesdarcyER-valid.pkl', 'rb') as input:
# BayesValid_ER = joblib.load(input)
sys.exit()
-
class Object(object):
pass
@@ -380,8 +385,8 @@ if __name__ == "__main__":
'logBME_ffpm-stokesdarcyER-valid.csv')
BayesValid_PNM.logBME = np.loadtxt(result_folder+'/outputs_ffpm-stokespnm/'+
'logBME_ffpm-stokespnm-valid.csv')
- BayesValid_PNM_NA.logBME = np.loadtxt(result_folder+'/outputs_ffpm-stokespnmNA/'+
- 'logBME_ffpm-stokespnmNA-valid.csv')
+ # BayesValid_PNM_NA.logBME = np.loadtxt(result_folder+'/outputs_ffpm-stokespnmNA/'+
+ # 'logBME_ffpm-stokespnmNA-valid.csv')
# BayesValid_ER.KLD = np.loadtxt('logBME_ER.csv')
# BayesValid_ER.infEntropy = np.loadtxt('logBME_ER.csv')
#BayesValid_PNM.logBME = np.loadtxt('logBME_PNM.csv')
@@ -393,7 +398,7 @@ if __name__ == "__main__":
BayesValid_PNM.BME = np.exp(BayesValid_PNM.logBME)
- BayesValid_PNM_NA.BME = np.exp(BayesValid_PNM_NA.logBME)
+ # BayesValid_PNM_NA.BME = np.exp(BayesValid_PNM_NA.logBME)
#==========================================================================
#============= Bayes Factor Computation (only models) =================
@@ -407,7 +412,7 @@ if __name__ == "__main__":
BME_Dict['Classical IC'] = BayesValid_BJ.BME
BME_Dict['Generalized IC'] = BayesValid_ER.BME
BME_Dict['Pore Network'] = BayesValid_PNM.BME
- BME_Dict['Pore Network SA'] = BayesValid_PNM_NA.BME
+ # BME_Dict['Pore Network SA'] = BayesValid_PNM_NA.BME
# kdePlot only models
kdePlot_BayesFactor(BME_Dict, plotName='Valid')
@@ -417,12 +422,14 @@ if __name__ == "__main__":
#==========================================================================
# ---------------- Validation ----------------
# Normalize the BME (Asumption: Model Prior weights are equal for models)
- ValidBME = np.vstack((BayesValid_BJ.logBME, BayesValid_ER.logBME, BayesValid_PNM.logBME, BayesValid_PNM_NA.logBME))
+ # ValidBME = np.vstack((BayesValid_BJ.logBME, BayesValid_ER.logBME, BayesValid_PNM.logBME, BayesValid_PNM_NA.logBME))
+ ValidBME = np.vstack((BayesValid_BJ.logBME, BayesValid_ER.logBME, BayesValid_PNM.logBME))
print(np.quantile(ValidBME, 0.25, axis=1)-np.quantile(ValidBME, 0.5, axis=1))
print(np.quantile(ValidBME, 0.5, axis=1))
print(np.quantile(ValidBME, 0.75, axis=1)-np.quantile(ValidBME, 0.5, axis=1))
- ValidBME = np.vstack((BayesValid_BJ.BME, BayesValid_ER.BME, BayesValid_PNM.BME, BayesValid_PNM_NA.BME))
+ # ValidBME = np.vstack((BayesValid_BJ.BME, BayesValid_ER.BME, BayesValid_PNM.BME, BayesValid_PNM_NA.BME))
+ ValidBME = np.vstack((BayesValid_BJ.BME, BayesValid_ER.BME, BayesValid_PNM.BME))
ValidModelWeightsBME = np.divide(ValidBME, np.nansum(ValidBME, axis=0))
print(np.quantile(ValidModelWeightsBME, 0.25, axis=1)-np.quantile(ValidModelWeightsBME, 0.5, axis=1))
diff --git a/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py b/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
index 4ab3ec3004ea2a7e50ca0b900d557da9593cf50a..9242132b7d014ebf63cb85f1c1604d3360affba9 100755
--- a/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
+++ b/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
@@ -59,7 +59,7 @@ def check_ranges(samples, BayesDF):
index.append(i)
return index
-def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
+def run(params, errorPerc=0.05, couplingcond='BJ', inletLoc='top', PCEEDMethod='normal'):
print("\n" + "="*75)
print("Stochastic calibration and validation of {0}.".format('ffpm-stokesdarcy'+couplingcond))
@@ -73,21 +73,23 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
Model = PyLinkForwardModel()
Model.nrCPUs = 4 #8
Model.Type = 'PyLink'
- Model.Name = 'ffpm-stokesdarcy'+couplingcond
- Model.InputFile = 'stokesdarcy_'+couplingcond+'_params.input'
- Model.InputTemplate = 'stokesdarcy_'+couplingcond+'_params.tpl.input'
+ Model.Name = 'ffpm-stokesdarcy{}_{}Open'.format(couplingcond,inletLoc)
+ Model.InputFile = 'stokesdarcy_{}_{}Open_params.input'.format(couplingcond,inletLoc)
+ Model.InputTemplate = 'stokesdarcy_{}_{}Open_params.tpl.input'.format(couplingcond,inletLoc)
+
+ Model.Command = "model_exe.py stokesdarcy {}".format(Model.InputFile)
- Model.Command = "model_exe.py stokesdarcy stokesdarcy_%s_params.input"%couplingcond
Model.ExecutionPath = os.getcwd()
Model.Output.Parser = 'read_ffpm'
Model.Output.Names = ['velocity [m/s]', 'p']
- Model.Output.FileNames = ["ffpm_stokesdarcy_velocity_final.csv","ffpm_stokesdarcy_p_final.csv"]
+ Model.Output.FileNames = ["ffpm_stokesdarcy_velocity_final.csv",
+ "ffpm_stokesdarcy_p_final.csv"]
# For Bayesian inversion
- Model.MeasurementFile = 'data/stokesDataCalib.csv'
+ Model.MeasurementFile = 'data/stokesDataCalib_{}Open.csv'.format(inletLoc)
# Include the validation observation data
- Model.MeasurementFileValid = 'data/stokesDataValid.csv'
+ Model.MeasurementFileValid = 'data/stokesDataValid_{}Open.csv'.format(inletLoc)
#=====================================================
#========= PROBABILISTIC INPUT MODEL ==============
@@ -162,8 +164,8 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
# The degree with the lowest Leave-One-Out cross-validation (LOO)
# error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
- MetaModelOpts.MinPceDegree = 5 # default = 1
- MetaModelOpts.MaxPceDegree = 5 #7
+ MetaModelOpts.MinPceDegree = 7 # default = 1
+ MetaModelOpts.MaxPceDegree = 7 #7
# q-quasi-norm 0<q<1 (default=1)
MetaModelOpts.q = 1.0 #np.linspace(0.3,0.6,3)
@@ -185,10 +187,11 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
# 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'
+ MetaModelOpts.ExpDesign.SamplingMethod = 'sobol'
# Provide the experimental design object with a hdf5 file
- # MetaModelOpts.ExpDesign.hdf5File = "ExpDesign_ffpm-stokesdarcy{}.hdf5".format(couplingcond)
+ # MetaModelOpts.ExpDesign.hdf5File = "ExpDesign_ffpm-stokesdarcy{}_{}Open.hdf5".format(couplingcond,inletLoc)
+
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
@@ -196,7 +199,7 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-16
# Defining the measurement error, if it's known a priori
- calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib.csv'))**2
+ calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib_{}Open.csv'.format(inletLoc)))**2
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
DiscrepancyOpts.Parameters = calibMeasuError
@@ -233,7 +236,7 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
# For colculation of validation error for SeqDesign
import h5py
- HDF5File = "ExpDesign_ffpm-stokesdarcy"+couplingcond+"-testset_Calibration.hdf5"
+ HDF5File = "ExpDesign_{}-testset_Calibration.hdf5".format(Model.Name)
ValidSets = h5py.File("./data/ValidationSets/"+HDF5File, 'r+')
MetaModelOpts.validSamples = np.array(ValidSets["EDX/init_"])
validModelRuns = dict()
@@ -322,12 +325,12 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
BayesOptsCalib.PlotMAPPred = False
# ----- Define the discrepancy model -------
- calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib.csv'))**2
+ calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib_{}Open.csv'.format(inletLoc)))**2
calibNumErr = pd.read_csv('data/NumError_stokesdarcy_'+couplingcond+'_Calib.csv')**2
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = calibMeasuError.add(calibNumErr,fill_value=0)
+ DiscrepancyOpts.Parameters = calibMeasuError#.add(calibNumErr,fill_value=0)
BayesOptsCalib.Discrepancy = DiscrepancyOpts
@@ -367,7 +370,8 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
#================= Visualization ===================
#=====================================================
# Plot posterior predictive
- postPredictiveplot(PCEModel.ModelObj.Name, errorPerc, case='Calib')
+ postPredictiveplot(PCEModel.ModelObj.Name, errorPerc, case='Calib',
+ inletLoc=inletLoc, bins=20)
#=====================================================
#================== VALIDATION =====================
@@ -397,19 +401,17 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
print("-"*40 + "\n")
import copy
ValidModel = copy.deepcopy(Model)
- ValidModel.Name = 'ffpm-stokesdarcy'+couplingcond+'-valid'
+ ValidModel.Name = '{}-valid'.format(Model.Name)
ValidMetaModelOpts = copy.deepcopy(MetaModelOpts)
ValidMetaModelOpts.Inputs = ValidInputs
ValidMetaModelOpts.ExpDesign.Method = 'normal'
ValidMetaModelOpts.ExpDesign.InputObj = ValidInputs
ValidMetaModelOpts.ExpDesign.Y = None
ValidMetaModelOpts.ExpDesign.SamplingMethod = 'random' #'random' 'user'
- # ValidMetaModelOpts.ExpDesign.hdf5File = "ExpDesign_ffpm-stokesdarcy{}-valid.hdf5".format(couplingcond)#None
- # ValidMetaModelOpts.MinPceDegree = 9 # default = 1
- # ValidMetaModelOpts.MaxPceDegree = 9
+ ValidMetaModelOpts.ExpDesign.hdf5File = None#"ExpDesign_{}-valid.hdf5".format(Model.Name)#None
import h5py
- HDF5File = "ExpDesign_ffpm-stokesdarcy"+couplingcond+"-testset_Validation.hdf5"
+ HDF5File = "ExpDesign_{}-testset_Validation.hdf5".format(Model.Name)
ValidSets = h5py.File("./data/ValidationSets/"+HDF5File, 'r+')
indices = check_ranges(np.array(ValidSets["EDX/init_"]), BayesCalib.Posterior_df)
ValidMetaModelOpts.validSamples = np.array(ValidSets["EDX/init_"])[indices]
@@ -457,12 +459,11 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
ValidPostPCE.seqDesignDiagnosticPlots()
# Compute and print RMSE error
- # ValidPCEModel.RMSE = None
try:
ValidPostPCE.accuracyCheckMetaModel(Samples=ValidMetaModelOpts.validSamples,
validOutputsDict=validModelRuns)
except:
- ValidPostPCE.accuracyCheckMetaModel(nSamples=5)
+ ValidPCEModel.RMSE = None
#=====================================================
#========= Bayesian inference (Validation) =========
@@ -489,12 +490,12 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
# # ----- Define the discrepancy model -------
- validMeasuError = (errorPerc * pd.read_csv('data/stokesDataValid.csv'))**2
+ validMeasuError = (errorPerc * pd.read_csv('data/stokesDataValid_{}Open.csv'.format(inletLoc)))**2
validNumErr = pd.read_csv('data/NumError_stokesdarcy_'+couplingcond+'_Valid.csv')**2
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = validMeasuError.add(validNumErr,fill_value=0)
+ DiscrepancyOpts.Parameters = validMeasuError#.add(validNumErr,fill_value=0)
BayesOptsValid.Discrepancy = DiscrepancyOpts
@@ -535,7 +536,8 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
#================= Visualization ===================
#=====================================================
# Plot posterior predictive
- postPredictiveplot(ValidPCEModel.ModelObj.Name, errorPerc, case='Valid')
+ postPredictiveplot(ValidPCEModel.ModelObj.Name, errorPerc, case='Valid',
+ inletLoc=inletLoc, bins=20)
#=====================================================
#======== Moving folders to a new folder ===========
@@ -590,4 +592,4 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
#==================== Run main scripts for PA-B =======================
#==========================================================================
# PCEModel, BayesCalib, BayesValid = run(params,couplingcond='BJ', PCEEDMethod=PCEExpDesignMethod)
-# PCEModel = run(params,couplingcond='BJ', PCEEDMethod=PCEExpDesignMethod)
+# PCEModel = run(params,couplingcond='ER',errorPerc=0.05,inletLoc='left' ,PCEEDMethod=PCEExpDesignMethod)
diff --git a/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py b/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
index ab0ad3e4d9068d2a8590e3b7ddc0c4bbc571bc50..f3a114ed4618846434f77ec3e806a6296bb4b23a 100755
--- a/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
+++ b/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
@@ -58,7 +58,7 @@ def check_ranges(samples, BayesDF):
index.append(i)
return index
-def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
+def run(params, averaging=True,errorPerc=0.05, inletLoc='top', PCEEDMethod='normal'):
print("\n" + "="*75)
print("Stochastic calibration and validation of {0}.".format('ffpm-stokespnm'))
@@ -67,11 +67,11 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
nInitSamples,nTotalSamples,nBootstrapItr,perturbedData = params
if averaging:
- calib_data_csv_path = 'data/stokesDataCalib.csv'
- valid_data_csv_path = 'data/stokesDataValid.csv'
+ calib_data_csv_path = 'data/stokesDataCalib_{}Open.csv'.format(inletLoc)
+ valid_data_csv_path = 'data/stokesDataValid_{}Open.csv'.format(inletLoc)
else:
- calib_data_csv_path = 'data/stokesDataCalib_without_averaging.csv'
- valid_data_csv_path = 'data/stokesDataValid_without_averaging.csv'
+ calib_data_csv_path = 'data/stokesDataCalib_{}Open_without_averaging.csv'.format(inletLoc)
+ valid_data_csv_path = 'data/stokesDataValid_{}Open_without_averaging.csv'.format(inletLoc)
#=====================================================
#============ COMPUTATIONAL MODEL ================
@@ -79,11 +79,14 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
Model = PyLinkForwardModel()
# Model.nrCPUs = 15
Model.Type = 'PyLink'
- Model.Name = 'ffpm-stokespnm' if averaging else 'ffpm-stokespnmNA'
- Model.InputFile = 'stokespnm_params.input'
- Model.InputTemplate = 'stokespnm_params.tpl.input'
+ if averaging:
+ Model.Name = 'ffpm-stokespnm_{}Open'.format(inletLoc)
+ else:
+ Model.Name = 'ffpm-stokespnmNA_{}Open'.format(inletLoc)
+ Model.InputFile = 'stokespnm_{}Open_params.input'.format(inletLoc)
+ Model.InputTemplate = 'stokespnm_{}Open_params.tpl.input'.format(inletLoc)
- Model.Command = "model_exe.py stokespnm stokespnm_params.input --averaging {}".format(averaging)
+ Model.Command = "model_exe.py stokespnm {} --averaging {}".format(Model.InputFile,averaging)
Model.ExecutionPath = os.getcwd()
Model.Output.Parser = 'read_ffpm'
Model.Output.Names = ['velocity [m/s]', 'p']
@@ -150,11 +153,11 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# The degree with the lowest Leave-One-Out cross-validation (LOO)
# error (or the highest score=1-LOO)estimator is chosen as the final
# metamodel.
- MetaModelOpts.MinPceDegree = 1 # default = 1 10
- MetaModelOpts.MaxPceDegree = 7 #10
+ MetaModelOpts.MinPceDegree = 5 # default = 1 10
+ MetaModelOpts.MaxPceDegree = 5 #10
# q-quasi-norm 0<q<1 (default=1)
- MetaModelOpts.q = 0.5 #np.linspace(0.3,0.6,3)
+ MetaModelOpts.q = 1.0 #np.linspace(0.3,0.6,3)
# Print summary of the regression results
# MetaModelOpts.DisplayFlag = True
@@ -176,7 +179,7 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
MetaModelOpts.ExpDesign.SamplingMethod = 'latin_hypercube'
# Provide the experimental design object with a hdf5 file
- # MetaModelOpts.ExpDesign.hdf5File = "ExpDesign_{}.hdf5".format(Model.Name)
+ # MetaModelOpts.ExpDesign.hdf5File = "ExpDesign_{}_{}open.hdf5".format(Model.Name,boundaryType)
# Sequential experimental design (needed only for sequential ExpDesign)
MetaModelOpts.ExpDesign.NrofNewSample = 1
@@ -275,7 +278,7 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# Compute and print RMSE error
PostPCE.accuracyCheckMetaModel(Samples=MetaModelOpts.validSamples,
validOutputsDict=validModelRuns)
-
+ # return PCEModel
#=====================================================
#========= Bayesian inference (Calibration) ========
#=====================================================
@@ -287,6 +290,7 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# Select the inference method
BayesOptsCalib.SamplingMethod = 'MCMC'
BayesOptsCalib.MCMCnwalkers = 300 #500
+ # BayesOptsCalib.selectedIndices = [4,5,6,7]
# import emcee
# BayesOptsCalib.MCMCmoves = emcee.moves.StretchMove()
# [(emcee.moves.DEMove(), 0.8),(emcee.moves.DESnookerMove(), 0.2),]
@@ -303,12 +307,12 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
BayesOptsCalib.PlotMAPPred = False
# ----- Define the discrepancy model -------
- calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib.csv'))**2
+ calibMeasuError = (errorPerc * pd.read_csv('data/stokesDataCalib_{}Open.csv'.format(inletLoc)))**2
calibNumErr = pd.read_csv('data/NumError_stokespnm_Calib.csv')**2
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = calibMeasuError.add(calibNumErr,fill_value=0)
+ DiscrepancyOpts.Parameters = calibMeasuError#.add(calibNumErr,fill_value=0)
BayesOptsCalib.Discrepancy = DiscrepancyOpts
# -- (Option C) --
@@ -346,7 +350,8 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
#================= Visualization ===================
#=====================================================
# Plot posterior predictive
- postPredictiveplot(PCEModel.ModelObj.Name, errorPerc, averaging, case='Calib')
+ postPredictiveplot(PCEModel.ModelObj.Name, errorPerc, averaging,
+ inletLoc=inletLoc, case='Calib', bins=20)
#=====================================================
#================== VALIDATION =====================
@@ -381,8 +386,8 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
ValidMetaModelOpts.ExpDesign.Method = 'normal'
ValidMetaModelOpts.ExpDesign.InputObj = ValidInputs
ValidMetaModelOpts.ExpDesign.Y = None
- ValidMetaModelOpts.ExpDesign.SamplingMethod = 'random'
- # ValidMetaModelOpts.ExpDesign.hdf5File = "ExpDesign_ffpm-stokespnm-valid.hdf5"#None
+ ValidMetaModelOpts.ExpDesign.SamplingMethod = 'random' #'random'
+ # ValidMetaModelOpts.ExpDesign.hdf5File = "ExpDesign_{}-valid.hdf5".format(Model.Name)#None
import h5py
HDF5File = "ExpDesign_{}-testset_Validation.hdf5".format(Model.Name)
@@ -437,7 +442,8 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
ValidPostPCE.accuracyCheckMetaModel(Samples=ValidMetaModelOpts.validSamples,
validOutputsDict=validModelRuns)
except:
- ValidPostPCE.accuracyCheckMetaModel(nSamples=5)
+ ValidPCEModel.RMSE = None
+ # ValidPostPCE.accuracyCheckMetaModel(nSamples=5)
#=====================================================
#========= Bayesian inference (Validation) =========
@@ -464,12 +470,12 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# # ----- Define the discrepancy model -------
- validMeasuError = (errorPerc * pd.read_csv('data/stokesDataValid.csv'))**2
+ validMeasuError = (errorPerc * pd.read_csv('data/stokesDataValid_{}Open.csv'.format(inletLoc)))**2
validNumErr = pd.read_csv('data/NumError_stokespnm_Valid.csv')**2
DiscrepancyOpts = Discrepancy('')
DiscrepancyOpts.Type = 'Gaussian'
- DiscrepancyOpts.Parameters = validMeasuError.add(validNumErr,fill_value=0)
+ DiscrepancyOpts.Parameters = validMeasuError#.add(validNumErr,fill_value=0)
BayesOptsValid.Discrepancy = DiscrepancyOpts
@@ -503,12 +509,13 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# Save class objects
with open('PA_A_Bayes'+ValidModel.Name+'.pkl', 'wb') as output:
joblib.dump(BayesValid, output, 2)
-
+ np.savetxt('logBME_'+ValidModel.Name+'.csv',BayesValid.logBME)
#=====================================================
#================= Visualization ===================
#=====================================================
# Plot posterior predictive
- postPredictiveplot(ValidPCEModel.ModelObj.Name, errorPerc, averaging, case='Valid')
+ postPredictiveplot(ValidPCEModel.ModelObj.Name, errorPerc, averaging,
+ case='Valid',inletLoc=inletLoc)
#=====================================================
#======== Moving folders to a new folder ===========
@@ -522,10 +529,12 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
metaModel = PCEModel if case == 'Calib' else ValidPCEModel
ModelName = metaModel.ModelObj.Name
files = ['ExpDesign_'+ModelName+'.hdf5',
- # ModelName + '.zip',
'PCEModel_'+ModelName+'.pkl',
'PA_A_Bayes'+ModelName+'.pkl']
+ if case == 'Valid':
+ files.append('logBME_'+ValidModel.Name+'.csv')
+
for file in files:
shutil.move("./"+file, outdir+'/'+file)
@@ -551,7 +560,9 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
# nTotalSamples = 200 #100 Total No. of orig. Model runs for surrogate training
# nBootstrapItr = 1000 # No. of bootstraping iterations for Bayesian analysis
# BootstrapNoise = 0.005 # Noise amount for bootstraping in Bayesian analysis
-# perturbedData = np.loadtxt('./data/perturbedValidData.csv',delimiter=',')
+# perturbedData = np.loadtxt('./data/perturbedValidDataAvg.csv',delimiter=',')
+# # perturbedData = np.loadtxt('./data/perturbedValidData.csv',delimiter=',')
+
# params = (nInitSamples,nTotalSamples,nBootstrapItr,perturbedData)
@@ -559,4 +570,4 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
#==================== Run main scripts for PA-B =======================
#==========================================================================
# PCEModel, BayesCalib, BayesValid = run(params,averaging=False, PCEEDMethod=PCEExpDesignMethod)
-# PCEModel = run(params,averaging=False, PCEEDMethod=PCEExpDesignMethod)
+# PCEModel = run(params,averaging=False, inletLoc='left', PCEEDMethod=PCEExpDesignMethod)
diff --git a/BayesValidRox/tests/PA-A/model_exe.py b/BayesValidRox/tests/PA-A/model_exe.py
index dee7f8d4d8e0e0ccc797859ebd8e269833daa2d8..9f9f6a8151d5bbfcce9ce07a4817e1b2eea71114 100755
--- a/BayesValidRox/tests/PA-A/model_exe.py
+++ b/BayesValidRox/tests/PA-A/model_exe.py
@@ -44,8 +44,8 @@ elif model == 'stokespnm':
model_dir = 'stokespnm'
pvdfile1 = 'regular_2D_stokespnm_250_micron_stokes'
pvdfile2 = 'regular_2D_stokespnm_250_micron_pnm'
- darcyfile = 'regular_2D_stokespnm_250_micron_pnm-00001.vtp'
- stokesfile = 'regular_2D_stokespnm_250_micron_stokes-00001.vtu'
+ darcyfile = 'regular_2D_stokespnm_250_micron_pnm-00000.vtp'
+ stokesfile = 'regular_2D_stokespnm_250_micron_stokes-00000.vtu'
velocity_points = "../models/velocity_points_pnm.csv"
pressure_points = "../models/pressure_points_pnm.csv"
@@ -72,7 +72,6 @@ if model == 'stokespnm':
-f "+pvdfile1+".pvd "+pvdfile2+".pvd\
-p "
-
Process2 = os.system(Command + velocity_points)
if Process2 != 0:
print('\nMessage 2:')
diff --git a/BayesValidRox/tests/PA-A/models/stokespnm/stokespnm_params.input b/BayesValidRox/tests/PA-A/models/stokespnm/stokespnm_params.input
deleted file mode 100644
index 364d36f9f5afd5af4bff71ffc354b6f9cf043169..0000000000000000000000000000000000000000
--- a/BayesValidRox/tests/PA-A/models/stokespnm/stokespnm_params.input
+++ /dev/null
@@ -1,62 +0,0 @@
-[Stokes.Grid]
-LowerLeft = -0.000125000 0.005125 # inlet length: 0.25e-2 m
-UpperRight = 0.010125 0.006125 # outlet length: 0.25e-2 m, heigth channel: 1e-3 m
-CellsPerThroat = 32
-Cells1 = 128
-FixedCellsBetweenThroats = 32
-Refinement = 0
-
-[PNM.Grid]
-LowerLeft = 0 125e-6
-UpperRight = 0.010 0.005125 # 20 x 10 block, block size 250e-6 m
-NumPores = 21 11
-PoreInscribedRadius = 125e-6
-ThroatInscribedRadius = 125e-6
-PriorityList = 2 3 0 1
-BoundaryFaceMarker = 1 1 1 2
-DeletionProbability = 0 0 1 1
-RemoveThroatsOnBoundary = 3
-ThroatCrossSectionShape = TwoPlates
-PoreGeometry = Cube
-
-[Problem]
-Name = regular_2D_stokespnm_250_micron
-EnableGravity = false
-
-[PNM.Problem]
-Name = pnm
-TransmissibilityThroat = 5.26480244078e-06
-TransmissibilityHalfPore = 1.7642875191945115e-05
-
-[Stokes.Problem]
-Name = stokes
-EnableInertiaTerms = false
-VyMaxTop = -1e-3
-TopInletLeft = 0.0015
-TopInletRight = 0.0045
-OutflowLength = 0.5e-3
-Beta = 35237.668793254459160 # 35238.834999009282910 with symm. grad.
-
-[Component]
-LiquidDensity = 1e3
-LiquidKinematicViscosity = 1e-6
-
-[Newton]
-MaxSteps = 20
-TargetSteps = 4
-MaxRelativeShift = 1e-10
-
-[Vtk]
-AddVelocity = 1
-WriteFaceData = false
-
-[FluxOverSurface]
-Verbose = false
-HorizontalPlaneY = 0.005
-
-[Assembly]
-NumericDifference.BaseEpsilon = 1e-4
-
-[FreeFlow]
-EnableUnsymmetrizedVelocityGradient = true
-EnableUnsymmetrizedVelocityGradientForBeaversJoseph = true
diff --git a/BayesValidRox/tests/PA-A/post_plot_PAPER.py b/BayesValidRox/tests/PA-A/post_plot_PAPER.py
old mode 100644
new mode 100755
index 5276b25584de76cfb60b6b32a456d742698ded11..86ba88d73b431b8503dd0be40959082c5c2c2aff
--- a/BayesValidRox/tests/PA-A/post_plot_PAPER.py
+++ b/BayesValidRox/tests/PA-A/post_plot_PAPER.py
@@ -52,16 +52,16 @@ def upper_rugplot(data, height=.05, ax=None, **kwargs):
lc = LineCollection(segs, transform=ax.get_xaxis_transform(), **kwargs)
ax.add_collection(lc)
-def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib', bins='auto'):
- result_folder = './Results_22_06_2021/outputs_{}/'.format(modelName.split('-v')[0])
- # result_folder = './'
+def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib',inletLoc='top', bins='auto'):
+ # result_folder = './outputs_{}/'.format(modelName.split('-v')[0])
+ result_folder = './'
directory = result_folder+'Outputs_Bayes_'+modelName+'_'+case
OutputDir = ('postPred_'+modelName+'_'+case)
if not os.path.exists(OutputDir): os.makedirs(OutputDir)
- #ave = '' if averaging else '_without_averaging'
- ave = '' #'_PNM_averaging' if 'ffpm-stokespnm' in modelName else ''
- data = pd.read_csv('./data/stokesData'+case+ave+'.csv')
+ ave = '' if averaging else '_without_averaging'
+ #ave = '' #'_PNM_averaging' if 'ffpm-stokespnm' in modelName else ''
+ data = pd.read_csv('data/stokesData{}_{}Open{}.csv'.format(case,inletLoc,ave))
# Load Post pred file
f = h5py.File(directory+'/'+"postPredictive.hdf5", 'r+')
@@ -104,7 +104,7 @@ def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib', bins=
# Posterior predictive
postPred = np.array(f["EDY/"+OutputName])[:,idx]
- sns.histplot(postPred[postPred>0], ax=ax, bins=bins,
+ sns.histplot(postPred, ax=ax, bins=bins, #[postPred>0]
color='orange',stat="count") #normalizes counts so that the sum of the bar heights is 1
# Reference data from the pore-scale simulation
@@ -155,6 +155,6 @@ def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib', bins=
fig.savefig('./'+OutputDir+'/'+plotname+'.pdf', bbox_inches='tight')
plt.close()
-modelName = 'ffpm-stokesdarcyBJ' #stokespnm stokesdarcyER stokesdarcyBJ
-postPredictiveplot(modelName, errorPrec=0.05, case='Calib')#, bins=20)
-postPredictiveplot(modelName+'-valid', errorPrec=0.05, case='Valid')#,bins=20)
+modelName = 'ffpm-stokesdarcyER_left' #stokespnm stokesdarcyER stokesdarcyBJ
+postPredictiveplot(modelName, errorPrec=0.05, averaging=True, case='Calib',inletLoc='left', bins=20)
+# postPredictiveplot(modelName+'-valid', errorPrec=0.05, case='Valid',bins=20)