You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

189 lines
4.9 KiB

import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle
import argparse
import pickle
import yaml
def is_ipython():
''' Check if script is run in IPython.
Returns:
bool: True if IPython, else False '''
try:
get_ipython()
ipy = True
except NameError:
ipy = False
return ipy
def load_data(file):
''' Load numpy data from file.
Returns
dict: data dictionary
'''
dat = np.load(file)
return dat
def plot_parameters(dat, input_file, deparameterize=False, ref=None):
''' Plot the parameters in separate subplots with uncertainties.
Args:
dat (dict): data dictionary
deparameterize (bool): flag indicating if parameters should be
deparameterized via 2**theta
ref: reference value to be plotted with parameters
'''
if is_ipython():
plt.ion()
idx_a = input_file.find('/')
idx_b = input_file[idx_a+1::].find('/')
name_file = input_file[idx_a+1:idx_b+idx_a+1]
inputfile_path = 'results/' + name_file + '/input.yaml'
with open(inputfile_path) as file:
inputfile = yaml.full_load(file)
#true_values = {
# 3: 3400,
# 4: 4200,
# 5: 11000,
# 6: 7800,
# 2: 100
# }
true_values = {
3: 4800,
4: 7020,
5: 11520,
6: 11520,
2: 75
}
true_values_c = {
3: 0.0008,
4: 0.00034,
5: 0.00034,
6: 0.00034,
2: 100
}
true_values_rp = {
3: 10,
4: 60,
5: 220,
6: 160,
2: 100
}
current_val = []
labels = []
ids = []
for bnd_c in inputfile['estimation']['boundary_conditions']:
if 'windkessel' in bnd_c['type']:
for bnd_set in inputfile['boundary_conditions']:
if bnd_c['id'] == bnd_set['id']:
ids.append(bnd_c['id'])
current_val.append(bnd_set['parameters']['R_d'])
labels.append('$R_' + str(bnd_c['id']))
elif 'dirichlet' in bnd_c['type']:
current_val.append(inputfile['boundary_conditions'][0]['parameters']['U'])
ids.append(bnd_c['id'])
labels.append('$U')
dim = dat['theta'].shape[-1]
fig1, axes = plt.subplots(1,1,figsize=(8,6))
axes.set_ylabel(r'$\theta$',fontsize=18)
t = dat['times']
theta = dat['theta']
P = dat['P_theta']
col = cycle(['C0', 'C1', 'C2', 'C3','C4'])
ls = cycle(['-', '-', '--', '--', ':', ':', '-.', '-.'])
#legends = cycle(['$R_3$','$R_4$','$R_5$','$R_6$','$U$'])
legends = cycle(labels)
col_ = next(col)
ls_ = next(ls)
legends_=next(legends)
if dim == 1:
theta = theta.reshape((-1, 1))
P = P.reshape((-1, 1, 1))
for i in range(dim):
true_level = np.log(true_values[ids[i]]/current_val[i])/np.log(2)
rec_value = np.round(2**theta[-1, i]*current_val[i],2)
cur_key = ids[i]
axes.plot(t, theta[:, i] + 1.5*i, '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$')
axes.fill_between(t, theta[:, i] + 1.5*i - np.sqrt(P[:, i, i]),
theta[:, i] + 1.5*i + np.sqrt(P[:, i, i]), alpha=0.3,
color=col_)
axes.plot(t,1.5*i + t*0 + true_level , color=col_,ls='--')
col_ = next(col)
legends_=next(legends)
axes.legend(fontsize=14,loc='lower right')
axes.set_xlim([-0.01,0.81])
axes.set_xlabel(r'time (s)',fontsize=18)
# print('theta_peak: \t {}'.format(theta[round(len(theta)/2), :]))
print('Final value theta: \t {}'.format(theta[-1, :]))
print('Deparameterized: 2^theta_end: \t {}'.format(2**theta[-1, :]))
print('Real values: \t {}'.format(true_values))
print('Recon values: \t {a}:{b} '.format(a=ids[:],b=np.round(2**theta[-1, :]*current_val,2)))
plt.savefig('windk_res')
if not is_ipython():
plt.show()
def get_parser():
parser = argparse.ArgumentParser(
description='''
Plot the time evolution of the ROUKF estimated parameters.
To execute in IPython::
%run plot_roukf_parameters.py [-d] [-r N [N \
...]] file
''',
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('file', type=str, help='path to ROUKF stats file')
parser.add_argument('-d', '--deparameterize', action='store_true',
help='deparameterize the parameters by 2**theta')
parser.add_argument('-r', '--ref', metavar='N', nargs='+', default=None,
type=float, help='Reference values for parameters')
return parser
if __name__ == '__main__':
args = get_parser().parse_args()
dat = load_data(args.file)
plot_parameters(dat, args.file,deparameterize=args.deparameterize, ref=args.ref)