P283370/NuMRI
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

asd

Browse Source
This commit is contained in:
J.E. Garay Labra 2021-10-21 15:50:32 +02:00
parent 4fae3807fd
commit 06e9bc2468
3 changed files with 61 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
kalman/graphics/figure3.py
View File

@ -50,10 +50,14 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
if is_ipython(): if is_ipython():
plt.ion() plt.ion()
resultsname = 'results_slices_15dB/'
idx_a = input_file.find('/') idx_a = input_file.find('/')
idx_b = input_file[idx_a+1::].find('/') idx_b = input_file[idx_a+1::].find('/')
name_file = input_file[idx_a+1:idx_b+idx_a+1] name_file = input_file[idx_a+1:idx_b+idx_a+1]
inputfile_path = 'results/' + name_file + '/input.yaml' inputfile_path = resultsname + name_file + '/input.yaml'
with open(inputfile_path) as file: with open(inputfile_path) as file:
@ -85,6 +89,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
else: else:
RC_flag = False RC_flag = False
line_split = 1.5 line_split = 1.5
rec_values = {}
current_val = [] current_val = []
current_val_C = [] current_val_C = []
ids_type = [] ids_type = []
@ -153,7 +159,7 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
std_down = 2**(-np.sqrt(P[:, idx, idx]))*curve std_down = 2**(-np.sqrt(P[:, idx, idx]))*curve
std_up = 2**np.sqrt(P[:, idx, idx])*curve std_up = 2**np.sqrt(P[:, idx, idx])*curve
dash_curve = true_values[ids[i]] + t*0 dash_curve = true_values[ids[i]] + t*0
rec_values[cur_key] = rec_value
@ -170,7 +176,7 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
axes3.set_box_aspect(1/4) axes3.set_box_aspect(1/4)
plt.xticks(fontsize=28) plt.xticks(fontsize=28)
plt.yticks(fontsize=28) plt.yticks(fontsize=28)
plt.savefig('results/' + name_file + '/U.png') plt.savefig(resultsname + name_file + '/U.png')
else: else:
axes1.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 3) axes1.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 3)
@ -223,7 +229,35 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
axes2.set_xlabel(r'$t (s)$',fontsize=36) axes2.set_xlabel(r'$t (s)$',fontsize=36)
fig2.savefig('C.png') fig2.savefig('C.png')
fig1.savefig('results/' + name_file + '/Rd.png') fig1.savefig(resultsname + name_file + '/Rd.png')
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 {}'.format(rec_values))
# saving reconstructed values
np.save(resultsname + name_file + '/recon_values.npy',rec_values)
print('----- paper ----')
rec_values_new = {}
epsilon = []
for rk in rec_values.keys():
eps = 100*rec_values[rk]/true_values[rk]-100
epsilon.append(eps)
if rk == 2:
rec_values_new[rk] = np.round(rec_values[rk],1)
else:
rec_values_new[rk] = np.round(rec_values[rk]/1000,2)
print('Recon. values: \t {}'.format(rec_values_new))
print('epsilon =' , np.round(np.mean(epsilon),2))
if not is_ipython(): if not is_ipython():
plt.show() plt.show()

18
kalman/graphics/figureRd.py
View File

@ -53,7 +53,7 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
idx_a = input_file.find('/') idx_a = input_file.find('/')
idx_b = input_file[idx_a+1::].find('/') idx_b = input_file[idx_a+1::].find('/')
name_file = input_file[idx_a+1:idx_b+idx_a+1] name_file = input_file[idx_a+1:idx_b+idx_a+1]
inputfile_path = 'results/' + name_file + '/input.yaml' inputfile_path = 'results_slices_15dB/' + name_file + '/input.yaml'
with open(inputfile_path) as file: with open(inputfile_path) as file:
inputfile = yaml.full_load(file) inputfile = yaml.full_load(file)
@ -107,9 +107,9 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
theta = dat['theta'] theta = dat['theta']
P = dat['P_theta'] P = dat['P_theta']
#col = cycle(['C0', 'C1', 'C2', 'C3','C4'])
col = cycle(['orangered', 'dodgerblue', 'limegreen', 'C3','C4']) color_list = ['tomato', 'springgreen' , '#2CBDFE']
col = cycle(color_list)
legends = cycle(labels) legends = cycle(labels)
col_ = next(col) col_ = next(col)
@ -132,7 +132,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
if ids_type[i] == 'dirichlet': if ids_type[i] == 'dirichlet':
axes3.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 5) #axes3.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 5)
axes3.plot(t, curve , '-', color=col_,label= legends_ + '$', linewidth = 5)
axes3.fill_between(t, std_down, std_up, alpha=0.3, color=col_) axes3.fill_between(t, std_down, std_up, alpha=0.3, color=col_)
legends_=next(legends) legends_=next(legends)
axes3.plot(t, dash_curve , color=col_,ls='--' , linewidth = 3) axes3.plot(t, dash_curve , color=col_,ls='--' , linewidth = 3)
@ -147,7 +148,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
#plt.savefig('U_' + name_file + '.png') #plt.savefig('U_' + name_file + '.png')
plt.close(fig3) plt.close(fig3)
else: else:
axes1.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 4) #axes1.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 4)
axes1.plot(t, curve , '-', color=col_,label= legends_ + '$', linewidth = 4)
axes1.fill_between(t, std_down, std_up, alpha=0.3, color=col_) axes1.fill_between(t, std_down, std_up, alpha=0.3, color=col_)
axes1.plot(t, dash_curve , color=col_,ls='--',linewidth = 3) axes1.plot(t, dash_curve , color=col_,ls='--',linewidth = 3)
legends_=next(legends) legends_=next(legends)
@ -161,8 +163,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
axes1.set_ylabel(r'$R_d$',fontsize=30) axes1.set_ylabel(r'$R_d$',fontsize=30)
axes1.legend(fontsize=36,loc='upper right') axes1.legend(fontsize=36,loc='upper right')
axes1.set_xlim([0,0.51]) axes1.set_xlim([0,0.8])
axes1.set_ylim([-1000,65000]) axes1.set_ylim([-1000,66000])
axes1.set_box_aspect(1/2) axes1.set_box_aspect(1/2)
plt.xticks(fontsize=28) plt.xticks(fontsize=28)
plt.yticks(fontsize=28) plt.yticks(fontsize=28)

38
kalman/input_files/aorta.yaml
View File

@ -8,7 +8,7 @@ fluid:
state_velocity: 'update' state_velocity: 'update'
io: io:
write_path: 'results/Pb_Hz2.3_SNR12V30' write_path: 'results/Pb_test'
restart: restart:
path: '' # './projects/nse_coa3d/results/test_restart2/' path: '' # './projects/nse_coa3d/results/test_restart2/'
time: 0 time: 0
@ -192,36 +192,24 @@ estimation:
measurements: measurements:
#-
# #mesh: '/home/yeye/NuMRI/kalman/meshes/coaortaH3_leo2.0.h5'
# mesh: '/home/yeye/Desktop/slices/slice_Hz5.7.h5'
# fe_degree: 0
# xdmf_file: 'measurements/slice_Hz5.7/Perturbation/Mg12V30/u_all.xdmf'
# file_root: 'measurements/slice_Hz5.7/Perturbation/Mg12V30/u{i}.h5'
# #xdmf_file: 'measurements/slice_Hz5.7/u_all.xdmf'
# #file_root: 'measurements/slice_Hz5.7/u{i}.h5'
# indices: 0 # indices of checkpoints to be processed. 0 == all
# velocity_direction: [0,0,1]
# #noise_stddev: 0 # standard deviation of Gaussian noise
# #noise_stddev: 22.39 # SNR 12 slice 5.7
# #noise_stddev: 15.76 # SNR 15 slice 5.7
# #noise_stddev: 8.75 # SNR 15 slice 2.3
# #noise_stddev: 12.15 # SNR 12 slice 2.3
# noise_stddev: 0.269
# VENC: 47
# module_meas_file_root: 'measurements/slice_Hz5.7/Perturbation/Mg12V30/module/M{i}.h5'
- -
#mesh: '/home/yeye/NuMRI/kalman/meshes/coaortaH3_leo2.0.h5'
mesh: '/home/yeye/Desktop/slices/slice_Hz2.3.h5' mesh: '/home/yeye/Desktop/slices/slice_Hz2.3.h5'
fe_degree: 0 fe_degree: 0
xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u_all.xdmf' xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg15V30/u_all.xdmf'
file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u{i}.h5' file_root: 'measurements/slice_Hz2.3/Perturbation/Mg15V30/u{i}.h5'
#xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg15V120/u_all.xdmf'
#file_root: 'measurements/slice_Hz2.3/Perturbation/Mg15V120/u{i}.h5'
#xdmf_file: 'measurements/slice_Hz2.3/u_all.xdmf'
#file_root: 'measurements/slice_Hz2.3/u{i}.h5'
indices: 0 # indices of checkpoints to be processed. 0 == all indices: 0 # indices of checkpoints to be processed. 0 == all
velocity_direction: [0,0,1] velocity_direction: [0,0,1]
noise_stddev: 0.25 #noise_stddev: 13.81
VENC: 26 noise_stddev: 'initial'
module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/module/M{i}.h5' #noise_stddev: 0.181
VENC: 28
module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg15V30/module/M{i}.h5'
roukf: roukf:
particles: 'simplex' # unique or simplex particles: 'simplex' # unique or simplex