From 79ca2180db19ffa82ecfec9ce16b9099327a1c89 Mon Sep 17 00:00:00 2001
From: farid <farid.mohammadi@iws.uni-stuttgart.de>
Date: Tue, 18 Jan 2022 15:21:19 +0100
Subject: [PATCH] [PA-A][util]add a post_plot with barchart for all models.
---
.../tests/PA-A/util/post_plot_PAPER_v2.py | 122 ++++++++++++++++++
1 file changed, 122 insertions(+)
create mode 100644 BayesValidRox/tests/PA-A/util/post_plot_PAPER_v2.py
diff --git a/BayesValidRox/tests/PA-A/util/post_plot_PAPER_v2.py b/BayesValidRox/tests/PA-A/util/post_plot_PAPER_v2.py
new file mode 100644
index 000000000..83c172bf9
--- /dev/null
+++ b/BayesValidRox/tests/PA-A/util/post_plot_PAPER_v2.py
@@ -0,0 +1,122 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb 18 08:28:13 2020
+
+@author: farid
+"""
+import numpy as np
+import os
+import sys
+import joblib
+import seaborn as sns
+import h5py
+import pandas as pd
+import matplotlib.pylab as plt
+
+# Add BayesValidRox path
+sys.path.insert(0, "./../../../../BayesValidRox/")
+plt.rcParams.update({'lines.markeredgewidth': 1})
+
+def postPredictiveplot(modelNames, result_folder= ".", case="Calib",
+ inclusion="squared", inletLoc="top", bins="auto"):
+
+ OutputDir = f"../{result_folder}/postPredPlots_{case}"
+ if not os.path.exists(OutputDir):
+ os.makedirs(OutputDir)
+
+ c_l = ['orange', 'royalblue', 'violet', 'forestgreen']
+ allMeanPred, allStdPred = {}, {}
+ point_IDs = {}
+ for outName in ['velocity [m/s]', 'p']:
+ allMeanPred[outName], allStdPred[outName] = {}, {}
+
+ for name in modelNames:
+ modelName = f'ffpm-{name}_{inclusion}_inclusion_{inletLoc}Inflow'
+
+ if case == 'calib':
+ directory = f"../{result_folder}/outputs_{modelName}/Outputs_Bayes_{modelName}_{case}"
+ else:
+ directory = f"../{result_folder}/outputs_{modelName}/Outputs_Bayes_{modelName}-valid_{case}"
+
+ # Load Post pred file
+ f = h5py.File(f"{directory}/postPredictive.hdf5", "r+")
+
+ # Load PCEModel
+ with open(f"../{result_folder}/outputs_{modelName}/PCEModel_{modelName}.pkl", "rb") as input:
+ PCEModel = joblib.load(input)
+
+ if case == "Calib":
+ data = pd.read_csv(f"../{PCEModel.ModelObj.MeasurementFile}")
+ else:
+ data = pd.read_csv(f"../{PCEModel.ModelObj.MeasurementFileValid}")
+
+ obs = data[outName][~np.isnan(data[outName])].values
+ # Generate Prior Predictive
+ # priorPred, std = PCEModel.eval_metamodel(nsamples=10000)
+
+ # Create a dataframe containing the predictions
+ # Add PointIds to the dict
+ point_IDs[outName] = np.array(f[f"x_values/{outName}"])
+
+ # Modify the name
+ if 'BJ' in name:
+ name = 'Classical IC'
+ elif 'ER' in name:
+ name = 'Generalized IC'
+ elif 'pnm' in name:
+ name = 'Pore Network'
+
+ # Add post predictions' mean and std to the dict
+ postPred = np.array(f[f"EDY/{outName}"])
+ allMeanPred[outName][name] = np.mean(postPred, axis=0)
+ allStdPred[outName][name] = np.std(postPred, axis=0)
+
+ # Add data to the dict
+ allMeanPred[outName]['Ref. data'] = obs
+ f.close()
+
+ # Plot barchart
+ for outName in ['velocity [m/s]', 'p']:
+ fig, ax = plt.subplots()
+ dict1 = {"Point ID": point_IDs[outName]}
+ dict2 = allMeanPred[outName]
+ right_order = ['Classical IC', 'Generalized IC', 'Pore Network',
+ 'Ref. data']
+ reordered_dict2 = {k: dict2[k] for k in right_order}
+ std = {k: allStdPred[outName][k] for k in right_order[:-1]}
+ df = pd.DataFrame({**dict1, **reordered_dict2})
+
+ width = 0.4 if outName == 'p' else 0.8
+
+ df.plot(x="Point ID", y=right_order, kind="bar", ax=ax, width=width,
+ yerr=std, capsize=5, rot=0, edgecolor="black", lw=1, color=c_l)
+
+ if outName != "p":
+ ax.set_ylabel(outName)
+ else:
+ ax.set_ylabel("pressure [bar]")
+
+ # save the current figure
+ plotname = outName if outName == "p" else "velocity"
+ fig.savefig(
+ f"./{OutputDir}/{plotname}_barchart.svg", bbox_inches="tight"
+ )
+ fig.savefig(
+ f"./{OutputDir}/{plotname}_barchart.pdf", bbox_inches="tight"
+ )
+ plt.close()
+
+ return allMeanPred, allStdPred
+
+
+modelNames = ['stokesdarcyBJ', 'stokesdarcyER', 'stokespnm']
+inletLoc = 'top' # top or left
+inclusion = 'squared' # squared or circular
+path = "Results_09_01_2022_topInflow"
+case = 'Valid'
+
+allMeanPred, allStdPred = postPredictiveplot(modelNames, inclusion=inclusion,
+ inletLoc=inletLoc,
+ result_folder=path,
+ case='Valid')
--
GitLab