[PA-A][util]add a post_plot with barchart for all models.

BayesValidRox/tests/PA-A/util/post_plot_PAPER_v2.py

+#!/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')