From 7b1e70e14dd05aace24cd57d537d7e93bc924842 Mon Sep 17 00:00:00 2001
From: farid <farid.mohammadi@iws.uni-stuttgart.de>
Date: Thu, 8 Jul 2021 09:14:33 +0200
Subject: [PATCH] [PA-A][util] added sobol indices plot script and modified
postDist_visualization.
---
BayesValidRox/tests/PA-A/post_plot.py | 8 +-
BayesValidRox/tests/PA-A/post_plot_PAPER.py | 160 ++++++++++++++++++
.../tests/PA-A/util/postDist_visualization.py | 64 +++++--
.../tests/PA-A/util/sobol_indices.py | 62 +++++++
4 files changed, 280 insertions(+), 14 deletions(-)
create mode 100644 BayesValidRox/tests/PA-A/post_plot_PAPER.py
create mode 100644 BayesValidRox/tests/PA-A/util/sobol_indices.py
diff --git a/BayesValidRox/tests/PA-A/post_plot.py b/BayesValidRox/tests/PA-A/post_plot.py
index 35b4b0c89..7ca0e7779 100755
--- a/BayesValidRox/tests/PA-A/post_plot.py
+++ b/BayesValidRox/tests/PA-A/post_plot.py
@@ -53,8 +53,8 @@ def upper_rugplot(data, height=.05, ax=None, **kwargs):
ax.add_collection(lc)
def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib', bins='auto'):
- #result_folder = './Results_20_05_2021/outputs_{}/'.format(modelName.split('-v')[0])
- result_folder = './'
+ result_folder = './Results_22_06_2021/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)
@@ -142,6 +142,6 @@ def postPredictiveplot(modelName, errorPrec, averaging=True, case='Calib', bins=
fig.savefig('./'+OutputDir+'/PointID_'+str(idx+1)+'_'+plotname+'.pdf', bbox_inches='tight')
plt.close()
-# modelName = 'ffpm-stokespnm' #stokespnm stokesdarcyER stokesdarcyBJ
-# postPredictiveplot(modelName, errorPrec=0.05, case='Calib', bins=75)
+modelName = 'ffpm-stokesdarcyER' #stokespnm stokesdarcyER stokesdarcyBJ
+postPredictiveplot(modelName, errorPrec=0.05, case='Calib', bins=75)
#postPredictiveplot(modelName+'-valid', errorPrec=0.05, case='Valid',bins=75)
diff --git a/BayesValidRox/tests/PA-A/post_plot_PAPER.py b/BayesValidRox/tests/PA-A/post_plot_PAPER.py
new file mode 100644
index 000000000..2f59ecfb0
--- /dev/null
+++ b/BayesValidRox/tests/PA-A/post_plot_PAPER.py
@@ -0,0 +1,160 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb 18 08:28:13 2020
+
+@author: farid
+"""
+import numpy as np
+import os, sys
+import seaborn as sns
+import h5py
+from scipy.stats import norm
+import pandas as pd
+
+try:
+ import cPickle as pickle
+except ModuleNotFoundError:
+ import pickle
+from matplotlib.patches import Patch
+import matplotlib.lines as mlines
+from matplotlib import ticker
+import matplotlib.pylab as plt
+# Add BayesValidRox path
+sys.path.insert(0,'./../../../../BayesValidRox/')
+plt.style.use('seaborn-deep') # talk, paper, poster
+#fsize = 50
+#params = {'legend.fontsize': fsize,
+ #'axes.labelsize': fsize,
+ #'axes.titlesize': fsize,
+ #'xtick.labelsize' :fsize,
+ #'ytick.labelsize': fsize,
+ #'grid.color': 'k',
+ #'grid.linestyle': ':',
+ #'grid.linewidth': 0.5,
+## 'mathtext.fontset' : 'stix',
+## 'mathtext.rm' : 'serif',
+ #'font.family' : 'serif',
+ #'font.serif' : "Times New Roman", # or "Times"
+ #'figure.figsize' : (32, 24),
+ ## 'savefig.dpi':1500,
+ #'text.usetex': True
+ #}
+#plt.rcParams.update(params)
+
+def upper_rugplot(data, height=.05, ax=None, **kwargs):
+ from matplotlib.collections import LineCollection
+ ax = ax or plt.gca()
+ kwargs.setdefault("linewidth", 1)
+ segs = np.stack((np.c_[data, data],
+ np.c_[np.ones_like(data), np.ones_like(data)-height]),
+ axis=-1)
+ 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 = './'
+ 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')
+
+ # Load Post pred file
+ f = h5py.File(directory+'/'+"postPredictive.hdf5", 'r+')
+
+ # Load PCEModel
+ with open(result_folder+'PCEModel_'+modelName+'.pkl', 'rb') as input:
+ PCEModel = pickle.load(input)
+
+ # Generate Prior Predictive
+ #priorPred, std = PCEModel.eval_metamodel(nsamples=10000)
+
+
+ for OutputName in ['velocity [m/s]', 'p']:
+ # x_coords = np.array(f["x_values/"+OutputName])
+
+ # Find pointIDs
+ csv_file = 'pressure_points.csv' if OutputName == 'p' else 'velocity_points.csv'
+ if case == 'Calib':
+ pointIDs = pd.read_csv('./models/'+csv_file).query("__vtkIsSelected__ \
+ == 'Calibration'")['vtkOriginalPointIds'].to_numpy()
+ else:
+ pointIDs = pd.read_csv('./models/'+csv_file).query("__vtkIsSelected__ \
+ == 'Validation'")['vtkOriginalPointIds'].to_numpy()
+ import matplotlib
+ gs00 = matplotlib.gridspec.GridSpec(3, 4)
+ fig = plt.figure()
+ cnt=0
+ for idx, x in enumerate(pointIDs):
+
+ if idx!=0 and idx%4==0: cnt += 1
+
+ if idx > 7:
+ ax = fig.add_subplot(gs00[cnt,idx%4+1:idx%4+2])
+ else:
+ ax = fig.add_subplot(gs00[cnt,idx%4])
+ # Prior predictive
+ #outputs = priorPred[OutputName][:,idx]
+ #sns.histplot(outputs[outputs>0], ax=ax,# bins=50,
+ #color='blue', alpha=0.4,stat="count")
+
+ # Posterior predictive
+ postPred = np.array(f["EDY/"+OutputName])[:,idx]
+ sns.histplot(postPred[postPred>0], ax=ax, bins=bins,
+ color='orange',stat="count") #normalizes counts so that the sum of the bar heights is 1
+
+ # Reference data from the pore-scale simulation
+ ax.axvline(x=data[OutputName][idx], linewidth=8, color='green')
+ #sns.histplot(np.random.normal(data[OutputName][idx], errorPrec*data[OutputName][idx], len(postPred[postPred>0])), ax=ax,bins=bins,
+ #color='green', alpha=0.5, stat="count")
+
+ # Print conidence interval of ExpDesign.Y (Trained area)
+ modelRuns = PCEModel.ExpDesign.Y[OutputName][:,idx]
+ upper_rugplot(modelRuns, ax=ax, alpha=0.75, color='grey')
+
+ # fmt = ticker.StrMethodFormatter("{x}")
+ # ax.xaxis.set_major_formatter(fmt)
+ # ax.yaxis.set_major_formatter(fmt)
+
+ legend_elements = [
+ Patch(facecolor='orange', edgecolor='orange',label='Posterior Pred.'),
+ Patch(facecolor='green', edgecolor='green',alpha=0.5, label='Ref. Data'),
+ mlines.Line2D([], [], marker='|', color='grey', alpha=0.75,
+ linestyle='None', markersize=15, markeredgewidth=1.5, label='Orig. Responses')]
+
+ font = {'family': 'serif',
+ 'weight': 'normal',
+ 'size': 28,
+ }
+
+ ax.legend(handles=legend_elements, fontsize=font['size']-12)
+
+
+ #x_label = 'Pressure [Pa]' if OutputName == 'p' else 'velocity [m/s]'
+ #ax.set_xlabel(x_label, fontdict=font)
+ ax.set_ylabel('Count', fontdict=font)
+ # ax.set_xscale('log')
+ ax.tick_params(axis='both', which='major', labelsize=font['size'])
+ plt.ticklabel_format(axis="x", style="sci", scilimits=(0,0))
+ plt.ticklabel_format(axis="y", style="sci", scilimits=(0,0))
+ ax.yaxis.get_offset_text().set_fontsize(font['size'])
+ ax.xaxis.get_offset_text().set_fontsize(font['size'])
+ plt.grid(True)
+
+ title = 'Point ID: '+str(x)
+ plt.title(title,fontdict=font)
+ fig.subplots_adjust(top=0.95)
+
+ plt.tight_layout(pad=0.2)
+ plotname = OutputName if OutputName == 'p' else 'velocity'
+ fig.savefig('./'+OutputDir+'/'+plotname+'.svg', bbox_inches='tight')
+ 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)
diff --git a/BayesValidRox/tests/PA-A/util/postDist_visualization.py b/BayesValidRox/tests/PA-A/util/postDist_visualization.py
index 0e6bcc617..72a407392 100755
--- a/BayesValidRox/tests/PA-A/util/postDist_visualization.py
+++ b/BayesValidRox/tests/PA-A/util/postDist_visualization.py
@@ -8,6 +8,7 @@ Created on Thu Aug 13 09:53:11 2020
import sys, os
import numpy as np
import pandas as pd
+import h5py
import corner
try:
import cPickle as pickle
@@ -17,15 +18,17 @@ except ModuleNotFoundError:
# Local
sys.path.insert(0,'./../../../../BayesValidRox/')
-modelName = 'stokesdarcyER' #stokespnm stokesdarcyER stokesdarcyBJ
+modelName = 'stokespnm' #stokespnm stokesdarcyER stokesdarcyBJ
# Load the pickle objects
-data_dir = '../Results_14_06_2021/outputs_ffpm-{}/'.format(modelName)
+data_dir = '../Results_22_06_2021/outputs_ffpm-{}/'.format(modelName)
+
+with open(data_dir+'PCEModel_ffpm-{}.pkl'.format(modelName), 'rb') as input:
+ PCEModel = pickle.load(input)
with open(data_dir+'PA_A_Bayesffpm-{}.pkl'.format(modelName), 'rb') as input:
Bayes = pickle.load(input)
-
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
@@ -53,7 +56,7 @@ plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
#=====================================================
#========== DEFINITION OF THE METAMODEL ============
#=====================================================
-PCEModel = Bayes.PCEModel
+# PCEModel = Bayes.PCEModel
# Update slicing index
# Extract the experimental design
EDX = PCEModel.ExpDesign.X
@@ -61,8 +64,8 @@ EDYDict = PCEModel.ExpDesign.Y
JDist = PCEModel.ExpDesign.JDist
X_train = PCEModel.ExpDesign.X
EDY = PCEModel.ExpDesign.Y
-
-
+Outputs = ['velocity [m/s]', 'p']
+labels = [PCEModel.Inputs.Marginals[i].Name for i in range(PCEModel.NofPa)]
#=====================================================
#================= Visualization ===================
#=====================================================
@@ -71,12 +74,53 @@ outDir = './posteriorPlot'
if not os.path.exists(outDir): os.makedirs(outDir)
# ----- Posterior plot ---------
-postSamples = Bayes.Posterior_df.to_numpy()
+# import emcee
+# path = '../Results_22_06_2021/outputs_ffpm-{}/Outputs_Bayes_ffpm-{}_Calib/emcee_sampler.h5'.format(modelName,modelName)
+# sampler = emcee.backends.HDFBackend(path)
+# try:
+# tau = sampler.get_autocorr_time(tol=0)
+# except emcee.autocorr.AutocorrError:
+# tau = 5
-# Plot with cornerplot
-labels = [Bayes.PCEModel.Inputs.Marginals[i].Name for i in range(PCEModel.NofPa)]
+# if all(np.isnan(tau)): tau = 5
+
+# burnin = int(2*np.nanmax(tau))
+# thin = int(0.5*np.nanmin(tau)) if int(0.5*np.nanmin(tau))!=0 else 1
+# postSamples = sampler.get_chain(discard=burnin, flat=True, thin=thin)
+
+# Remove bad chains
+# limit = {1: 0.0007, 2:2.5e-5, 3:4e-8} # pnm
+limit = {1: 0.00088, 2:5e-3, 3:4e-8} #{1: 0.0007, 2:0, 3:4e-8} ER
-figPosterior = corner.corner(postSamples, labels=labels,
+# postSamples = postSamples.reshape(-1, postSamples.shape[-1])
+postSamples = Bayes.Posterior_df.to_numpy()
+goodIndices = (postSamples[:,1]>limit[2]) #(postSamples[:,0]>limit[1]) & (postSamples[:,1]>limit[2])#&
+ # (postSamples[:,2]>limit[3])]
+goodpostSamples = postSamples#[goodIndices]
+
+# Save new Posterior_df in Bayes.pkl file
+# Bayes.Posterior_df = pd.DataFrame(goodpostSamples, columns=labels)
+# bayesDir = 'Outputs_Bayes_ffpm-{}_Calib/'.format(modelName)
+# with open(data_dir+'PA_A_Bayesffpm-{}.pkl'.format(modelName), 'wb') as output:
+# pickle.dump(Bayes, output, pickle.HIGHEST_PROTOCOL)
+
+# # Remove take good samples in postPredictive
+# # Load Post pred file
+# f = h5py.File(data_dir+bayesDir+'/'+"postPredictive.hdf5", 'r+')
+# for out in Outputs:
+# postPred = f["EDY/"+out]
+# newpostPred = np.array(postPred)[goodIndices]
+# del f["EDY/"+out]
+# dset = f.create_dataset("EDY/"+out, data=newpostPred)
+# f.close()
+
+if 'pnm' in modelName:
+ labels[1] = '$g_{t,ij}$'
+ labels[2] = '$g_{p,i}$'
+ labels[3] = '$\\beta_{pore}$'
+
+# Plot with cornerplot
+figPosterior = corner.corner(goodpostSamples, labels=labels,
# range=PCEModel.ExpDesign.BoundTuples,
# color='grey',
use_math_text = True,
diff --git a/BayesValidRox/tests/PA-A/util/sobol_indices.py b/BayesValidRox/tests/PA-A/util/sobol_indices.py
new file mode 100644
index 000000000..a40c1bf0f
--- /dev/null
+++ b/BayesValidRox/tests/PA-A/util/sobol_indices.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 24 17:28:39 2021
+
+@author: farid
+"""
+import sys
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+SIZE = 60
+plt.rc('figure', figsize = (24, 16))
+plt.rc('font', family='serif', serif='Arial')
+plt.rc('font', size=SIZE)
+plt.rc('axes', grid = True)
+plt.rc('text', usetex=True)
+plt.rc('axes', linewidth=3)
+plt.rc('axes', grid=True)
+plt.rc('grid', linestyle="-")
+plt.rc('axes', titlesize=SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=SIZE) # legend fontsize
+plt.rc('figure', titlesize=SIZE) # fontsize of the figure title
+
+result_dir = '../Results_22_06_2021/'
+modelName = 'stokespnm' #stokesdarcyER stokespnm
+
+
+fig = plt.figure()
+
+xlabel = 'Point ID'
+x_values = {'velocity [m_s]':[1,2,3,4,5,6,7,8,9,10], 'p': [1,2,3]}
+outNames = ['velocity [m_s]','p']
+
+
+for outIdx, Output in enumerate(outNames):
+
+ total_sobols = pd.read_csv(result_dir+'outputs_ffpm-{}/Outputs_PostProcessing_calib/totalsobol_{}.csv'.format(modelName,Output))
+ total_sobol = total_sobols.to_numpy()
+ parNames = list(total_sobols.keys())
+
+ if 'pnm' in modelName:
+ parNames[1] = '$g_{t,ij}$'
+ parNames[2] = '$g_{p,i}$'
+ parNames[3] = '$\\beta_{pore}$'
+
+ ax = fig.add_axes([0,0,1,1])
+ dict1 = {xlabel: x_values[Output]}
+ dict2 = {param: sobolIndices for param, sobolIndices in zip(parNames,total_sobol.T)}
+
+ df = pd.DataFrame({**dict1, **dict2})
+ df.plot(x=xlabel, y=parNames, kind="bar", ax=ax , rot=0, colormap='Dark2')
+ ax.set_ylabel('Total Sobol indices, $S^T$')
+
+ # save the current figure
+ fig.savefig('Sobol_indices_{}_{}.pdf'.format(Output,modelName), bbox_inches='tight')
+
+ # Destroy the current plot
+ plt.clf()
--
GitLab