diff --git a/kalman/graphics/figure3.py b/kalman/graphics/figure3.py index 4e2b411..8d80f15 100644 --- a/kalman/graphics/figure3.py +++ b/kalman/graphics/figure3.py @@ -50,10 +50,14 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None): if is_ipython(): plt.ion() + + + resultsname = 'results_slices_15dB/' + 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' + inputfile_path = resultsname + name_file + '/input.yaml' with open(inputfile_path) as file: @@ -85,6 +89,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None): else: RC_flag = False line_split = 1.5 + + rec_values = {} current_val = [] current_val_C = [] 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_up = 2**np.sqrt(P[:, idx, idx])*curve 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) plt.xticks(fontsize=28) plt.yticks(fontsize=28) - plt.savefig('results/' + name_file + '/U.png') + plt.savefig(resultsname + name_file + '/U.png') else: 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) 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(): plt.show() diff --git a/kalman/graphics/figureRd.py b/kalman/graphics/figureRd.py index 4d18b94..ac489d8 100644 --- a/kalman/graphics/figureRd.py +++ b/kalman/graphics/figureRd.py @@ -53,7 +53,7 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None): 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' + inputfile_path = 'results_slices_15dB/' + name_file + '/input.yaml' with open(inputfile_path) as file: inputfile = yaml.full_load(file) @@ -107,9 +107,9 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None): theta = dat['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) col_ = next(col) @@ -132,7 +132,8 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None): 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_) legends_=next(legends) 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.close(fig3) 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.plot(t, dash_curve , color=col_,ls='--',linewidth = 3) 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.legend(fontsize=36,loc='upper right') - axes1.set_xlim([0,0.51]) - axes1.set_ylim([-1000,65000]) + axes1.set_xlim([0,0.8]) + axes1.set_ylim([-1000,66000]) axes1.set_box_aspect(1/2) plt.xticks(fontsize=28) plt.yticks(fontsize=28) diff --git a/kalman/input_files/aorta.yaml b/kalman/input_files/aorta.yaml index 0b3b294..3f6dd7d 100755 --- a/kalman/input_files/aorta.yaml +++ b/kalman/input_files/aorta.yaml @@ -8,7 +8,7 @@ fluid: state_velocity: 'update' io: - write_path: 'results/Pb_Hz2.3_SNR12V30' + write_path: 'results/Pb_test' restart: path: '' # './projects/nse_coa3d/results/test_restart2/' time: 0 @@ -192,36 +192,24 @@ estimation: 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' fe_degree: 0 - xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u_all.xdmf' - file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u{i}.h5' + xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg15V30/u_all.xdmf' + 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 velocity_direction: [0,0,1] - noise_stddev: 0.25 - VENC: 26 - module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/module/M{i}.h5' + #noise_stddev: 13.81 + noise_stddev: 'initial' + #noise_stddev: 0.181 + VENC: 28 + module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg15V30/module/M{i}.h5' - roukf: particles: 'simplex' # unique or simplex