diff --git a/kalman/graphics/figureRd.py b/kalman/graphics/figureRd.py new file mode 100644 index 0000000..4d18b94 --- /dev/null +++ b/kalman/graphics/figureRd.py @@ -0,0 +1,205 @@ +import matplotlib.pyplot as plt +import numpy as np +from itertools import cycle +import argparse +import pickle +import yaml + + +#import matplotlib.font_manager +from matplotlib import rc +rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) +rc('text', usetex=True) + + + +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: 4800, + 4: 7200, + 5: 11520, + 6: 11520, + 2: 75 + } + + + dim = dat['theta'].shape[-1] + + current_val = [] + ids_type = [] + 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']) + ids_type.append('windkessel') + current_val.append(bnd_set['parameters']['R_d']) + labels.append('$R_' + str(bnd_c['id']-3)) + + + elif 'dirichlet' in bnd_c['type']: + current_val.append(inputfile['boundary_conditions'][0]['parameters']['U']) + ids.append(bnd_c['id']) + ids_type.append('dirichlet') + labels.append('$U') + + + if 'windkessel' in ids_type: + fig1, axes1 = plt.subplots(1,1,figsize=(12,7)) + if 'dirichlet' in ids_type: + fig3, axes3 = plt.subplots(1,1,figsize=(12,7)) + + + + + + t = dat['times'] + theta = dat['theta'] + P = dat['P_theta'] + + #col = cycle(['C0', 'C1', 'C2', 'C3','C4']) + col = cycle(['orangered', 'dodgerblue', 'limegreen', 'C3','C4']) + + legends = cycle(labels) + + col_ = next(col) + legends_=next(legends) + + if dim == 1: + theta = theta.reshape((-1, 1)) + P = P.reshape((-1, 1, 1)) + + idx = 0 + + for i in range(len(ids)): + + cur_key = ids[i] + rec_value = np.round(2**theta[-1, idx]*current_val[i],2) + curve = 2**theta[:, idx]*current_val[i] + 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 + + + if ids_type[i] == 'dirichlet': + axes3.plot(t, curve , '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', 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) + axes3.set_ylabel(r'$U$',fontsize=36) + axes3.legend(fontsize=36,loc='upper right') + axes3.set_xlim([0,0.45]) + axes3.set_ylim([8,180]) + axes3.set_xlabel(r'$t (s)$',fontsize=36) + axes3.set_box_aspect(1/2) + plt.xticks(fontsize=28) + plt.yticks(fontsize=28) + #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.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) + + + col_ = next(col) + idx +=1 + + + + + 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_box_aspect(1/2) + plt.xticks(fontsize=28) + plt.yticks(fontsize=28) + axes1.set_xlabel(r'$t (s)$',fontsize=36) + + path_paper = '/home/yeye/A_aliasing_kalman/latex/0_preprint/Figures/' + path_paper = '/home/yeye/Desktop/' + #fig1.savefig('Rd_'+ name_file +'.png') + fig1.savefig(path_paper + 'Rd_'+ name_file +'.png') + + + + + +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) diff --git a/kalman/graphics/figureU.py b/kalman/graphics/figureU.py new file mode 100644 index 0000000..25cc9d5 --- /dev/null +++ b/kalman/graphics/figureU.py @@ -0,0 +1,153 @@ +import matplotlib.pyplot as plt +import numpy as np +from itertools import cycle +import argparse +import pickle +import yaml + + +#import matplotlib.font_manager +from matplotlib import rc +rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) +rc('text', usetex=True) + + + + + +def plot_parameters(): + ''' 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 + ''' + + + + name_file = ['SNR12V120_Pf','SNR12V70_Pf','SNR12V30_Pf'] + #name_file = ['SNR12V120_Pf_MAG','SNR12V70_Pf_MAG','SNR12V30_Pf_MAG'] + #name_file = ['SNR12V120_Pb_MAG','SNR12V70_Pb_MAG','SNR12V30_Pb_MAG'] + + name_file = ['slice2.3_Pa'] + vencs = ['180','105','45'] + path0 = '/home/yeye/Desktop/kalman/results/' + + + fig, axes = plt.subplots(1,1,figsize=(12,7)) + col = cycle(['orangered', 'dodgerblue', 'limegreen', 'C3','C4']) + + true_values = { + 3: 4800, + 4: 7200, + 5: 11520, + 6: 11520, + 2: 75 + } + + + + for nn,name in enumerate(name_file): + + + + path1 = path0 + name + '/' + inputfile_path = path1 + 'input.yaml' + dat = np.load(path1 + 'theta_stats.npz') + + + with open(inputfile_path) as file: + inputfile = yaml.full_load(file) + + col_ = next(col) + + + dim = dat['theta'].shape[-1] + + current_val = [] + ids_type = [] + 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']) + ids_type.append('windkessel') + current_val.append(bnd_set['parameters']['R_d']) + + + + elif 'dirichlet' in bnd_c['type']: + current_val.append(inputfile['boundary_conditions'][0]['parameters']['U']) + ids.append(bnd_c['id']) + ids_type.append('dirichlet') + labels.append('$U') + + + + + t = dat['times'] + theta = dat['theta'] + P = dat['P_theta'] + + legends = cycle(labels) + + legends_=next(legends) + + if dim == 1: + theta = theta.reshape((-1, 1)) + P = P.reshape((-1, 1, 1)) + + idx = 0 + + for i in range(len(ids)): + + cur_key = ids[i] + rec_value = np.round(2**theta[-1, idx]*current_val[i],2) + curve = 2**theta[:, idx]*current_val[i] + 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 + + + if ids_type[i] == 'dirichlet': + axes.plot(t, curve , '-', color=col_,label= '$(venc \ '+ vencs[nn] + ' \ cm/s) \ U = ' + str(rec_value) + '/' + str(true_values[cur_key]) + '$', linewidth = 5) + #axes.plot(t, curve , '-', color=col_,label= legends_ + vencs[nn] + ' \ cm/s$', linewidth = 5) + axes.fill_between(t, std_down, std_up, alpha=0.3, color=col_) + legends_=next(legends) + axes.plot(t, dash_curve , color='black',ls='--' , linewidth = 3) + + + idx +=1 + + + + axes.set_ylabel(r'$U$',fontsize=36) + axes.legend(fontsize=30,loc='upper right') + axes.set_xlim([0,0.35]) + axes.set_ylim([10,160]) + axes.set_xlabel(r'$t (s)$',fontsize=36) + axes.set_box_aspect(1/2) + plt.xticks(fontsize=28) + plt.yticks(fontsize=28) + plt.savefig('U.png') + plt.show() + + + + #path_paper = '/home/yeye/A_aliasing_kalman/latex/0_preprint/Figures/' + #fig1.savefig('Rd_'+ name_file +'.png') + + + + + + +if __name__ == '__main__': + + plot_parameters() diff --git a/kalman/input_files/aorta.yaml b/kalman/input_files/aorta.yaml index 63a4893..0b3b294 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_Hz5.7' + write_path: 'results/Pb_Hz2.3_SNR12V30' restart: path: '' # './projects/nse_coa3d/results/test_restart2/' time: 0 @@ -146,7 +146,7 @@ fem: streamline_diffusion: enabled: False parameter: 'standard' # standard, shakib, codina, klr - length_scale: 'metric' # average, max, metric + length_scale: 'average' # average, max, metric parameter_element_constant: True Cinv: ~ monolithic: @@ -192,36 +192,39 @@ 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/Mg12V120/u_all.xdmf' - file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V120/u{i}.h5' - #xdmf_file: 'measurements/slice_Hz5.7/u_all.xdmf' - #file_root: 'measurements/slice_Hz5.7/u{i}.h5' + xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u_all.xdmf' + file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/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 - - - mesh: '/home/yeye/Desktop/slices/slice_Hz5.7.h5' - fe_degree: 0 - xdmf_file: 'measurements/slice_Hz5.7/Perturbation/Mg12V120/u_all.xdmf' - file_root: 'measurements/slice_Hz5.7/Perturbation/Mg12V120/u{i}.h5' - indices: 0 # indices of checkpoints to be processed. 0 == all - velocity_direction: [0,0,1] - noise_stddev: 22.39 # SNR 12 slice 5.7 + noise_stddev: 0.25 + VENC: 26 + module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/module/M{i}.h5' + roukf: particles: 'simplex' # unique or simplex observation_operator: 'postprocessing' #state or postprocessing reparameterize: True - MAG_functional: - enable: False - VENC: 61 - module_meas_file_root: 'measurements/aorta_zdir/Perturbation/Mg15V30/module/M{i}.h5' \ No newline at end of file + MAG_functional: True \ No newline at end of file