asd

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)

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()

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'
+            xdmf_file: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/u_all.xdmf'
-            file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V120/u{i}.h5'
+            file_root: 'measurements/slice_Hz2.3/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: 0.25
-            #noise_stddev: 22.39     # SNR 12 slice 5.7
+            VENC: 26
-            #noise_stddev: 15.76     # SNR 15 slice 5.7
+            module_meas_file_root: 'measurements/slice_Hz2.3/Perturbation/Mg12V30/module/M{i}.h5'
-            #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
             
     
     roukf:
         particles: 'simplex' # unique or simplex
         observation_operator: 'postprocessing'         #state or postprocessing
         reparameterize: True
-        MAG_functional:
+        MAG_functional: True
-            enable: False
-            VENC: 61
-            module_meas_file_root: 'measurements/aorta_zdir/Perturbation/Mg15V30/module/M{i}.h5'