diff --git a/linfa/tests/test_no_disc_TP15_error_estimation.py b/linfa/tests/test_no_disc_TP15_error_estimation.py index 6c65783..725fc5d 100644 --- a/linfa/tests/test_no_disc_TP15_error_estimation.py +++ b/linfa/tests/test_no_disc_TP15_error_estimation.py @@ -180,13 +180,12 @@ def log_prior(calib_inputs, transform): # Add gaussian log prior for first two parameters gauss_prior_res = l1 + l2 + l3 - # Add beta prior for third parameter if False: + # Add beta prior for third parameter sigma_prior = torch.distributions.beta.Beta(torch.tensor([1.0]), torch.tensor([0.5])) prior_res = sigma_prior.log_prob(phys_inputs[:,2]) - print(prior_res) - exit() else: + # Add uniform prior a = torch.tensor([0.0]) b = torch.tensor([0.15]) sigma_prior = torch.distributions.uniform.Uniform(low = a, high = b) @@ -195,11 +194,8 @@ def log_prior(calib_inputs, transform): if a <= siggy <= b: prior_res[loopA] = sigma_prior.log_prob(siggy) else: + # Large negative number prior_res[loopA] = 500 - - # print('noise param', phys_inputs[:,2]) - # print('prior results', prior_res) - # exit() res = gauss_prior_res + prior_res + adjust diff --git a/plot_mcmc_and_linfa.py b/plot_mcmc_and_linfa.py new file mode 100644 index 0000000..7d4d48d --- /dev/null +++ b/plot_mcmc_and_linfa.py @@ -0,0 +1,157 @@ +import os +import argparse +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.ticker as mtick +from matplotlib.ticker import ScalarFormatter, MaxNLocator +from scipy.stats import gaussian_kde + +# Set global plot settings +plt.rcParams['figure.figsize'] = (8, 6) +plt.rcParams['figure.dpi'] = 300 +plt.rcParams['axes.labelsize'] = 16 +plt.rcParams['xtick.labelsize'] = 15 +plt.rcParams['ytick.labelsize'] = 15 +plt.rcParams['legend.fontsize'] = 12 +plt.rcParams['lines.linewidth'] = 3 +plt.rcParams['lines.markersize'] = 16 +plt.rcParams['axes.labelweight'] = 'bold' +plt.rcParams['xtick.direction'] = 'in' +plt.rcParams['ytick.direction'] = 'in' +plt.rcParams['xtick.top'] = True +plt.rcParams['ytick.right'] = True +plt.rcParams['savefig.bbox'] = 'tight' + +def plot_params(param_data, LL_data, idx1, idx2, param_file_mcmc_1, param_file_mcmc_2, out_dir,out_info,fig_format='png', use_dark_mode=False): + + # Read data + param_data = np.loadtxt(param_data) # [5000, 3] + dent_data = np.loadtxt(LL_data) # [5000, ] + mcmc_1_data = np.loadtxt(param_file_mcmc_1) + mcmc_2_data = np.loadtxt(param_file_mcmc_2) + + # Combine MCMC samples + samples = np.vstack([mcmc_1_data, mcmc_2_data]) # Transpose to get shape (n, d) + + # Perform KDE + kde = gaussian_kde(samples) + + # Create a grid to evaluate KDE + x = np.linspace(samples[0].min(), samples[0].max(), 100) + y = np.linspace(samples[1].min(), samples[1].max(), 100) + X, Y = np.meshgrid(x, y) + positions = np.vstack([X.ravel(), Y.ravel()]) + Z = np.reshape(kde(positions), X.shape) + + # Plot + plt.figure(figsize=(8, 6)) + plt.contour(X, Y, Z) + plt.scatter(param_data[:,idx1], param_data[:,idx2], lw = 0, s = 7, marker = 'o', c = np.exp(dent_data)) + plt.plot(1000, -21.0E3, 'r*') + plt.colorbar() + plt.xlabel(r'$\theta_{K,'+str(idx1+1)+'}$') + plt.ylabel(r'$\theta_{K,'+str(idx2+1)+'}$') + plt.savefig(os.path.join(out_dir,'params_plot_' + out_info + '_'+str(idx1)+'_'+str(idx2)+'.'+fig_format)) + plt.close() + + plt.figure(figsize=(6, 6)) + plt.hist(param_data[:,idx1], color = 'blue', alpha = 0.25, label = 'LINFA', density = True) + plt.hist(mcmc_1_data, color = 'red', alpha = 0.25, label = 'MCMC', density = True) + plt.axvline(1000, color = 'r') + plt.xlabel(r'$\theta_{K,'+str(idx1+1)+'}$') + plt.legend() + plt.savefig(os.path.join(out_dir,'marginal_params_plot_' + out_info + '_'+str(idx1)+'_'+str(idx2)+'.'+fig_format)) + plt.close() + + plt.figure(figsize=(6, 6)) + plt.hist(param_data[:,idx2], color = 'blue', alpha = 0.25, label = 'LINFA', density = True) + plt.hist(mcmc_2_data, color = 'red', alpha = 0.25, label = 'MCMC', density = True) + plt.axvline(-21.0E3, color = 'r') + plt.xlabel(r'$\theta_{K,'+str(idx2+1)+'}$') + plt.legend() + plt.savefig(os.path.join(out_dir,'marginal_params_plot_' + out_info + '_'+str(idx2)+'_'+str(idx2)+'.'+fig_format)) + plt.close() + +# ========= +# MAIN CODE +# ========= +if __name__ == '__main__': + + # Init parser + parser = argparse.ArgumentParser(description='.') + + # folder name + parser.add_argument('-f', '--folder', + action=None, + const=None, + default='./', + type=str, + required=False, + help='Folder with experiment results', + metavar='', + dest='folder_name') + + # folder name + parser.add_argument('-n', '--name', + action=None, + const=None, + default='./', + type=str, + required=True, + help='Name of numerical experiment', + metavar='', + dest='exp_name') + + # iteration number = 1 + parser.add_argument('-i', '--iter', + action=None, + const=None, + default=1, + type=int, + choices=None, + required=True, + help='Iteration number', + metavar='', + dest='step_num') + + # plot format + parser.add_argument('-p', '--picformat', + action=None, + const=None, + default='png', + type=str, + choices=['png','pdf','jpg'], + required=False, + help='Output format for picture', + metavar='', + dest='img_format') + + # Enable dark mode for pictures + parser.add_argument('-d', '--dark', + action='store_true', + default=False, + required=False, + help='Generate pictures for dark background', + dest='use_dark_mode') + + # Parse Commandline Arguments + args = parser.parse_args() + + # Set file name/path + out_dir = os.path.join(args.folder_name,args.exp_name) + param_file = os.path.join(out_dir,args.exp_name + '_params_' + str(args.step_num)) + LL_file = os.path.join(out_dir,args.exp_name + '_logdensity_' + str(args.step_num)) + out_info = args.exp_name + '_' + str(args.step_num) + param_file_mcmc_1 = os.path.join(out_dir,'posterior_samples_1') + param_file_mcmc_2 = os.path.join(out_dir,'posterior_samples_2') + + # Plot 2D slice of posterior samples + if(os.path.isfile(param_file) and os.path.isfile(LL_file)): + tot_params = np.loadtxt(param_file).shape[1] # extract total number of parameters inferred + print('Plotting posterior samples...') + for loopA in range(tot_params): # loop over total number of parameters + for loopB in range(loopA+1, tot_params): # get next parameter + plot_params(param_file,LL_file,loopA,loopB,param_file_mcmc_1, param_file_mcmc_2, out_dir,out_info,fig_format=args.img_format,use_dark_mode=args.use_dark_mode) + else: + print('File with posterior samples not found: '+param_file) + print('File with log-density not found: '+LL_file) diff --git a/run_plot_res.sh b/run_plot_res.sh index b3674db..d39d0a7 100644 --- a/run_plot_res.sh +++ b/run_plot_res.sh @@ -1,7 +1,7 @@ -python3 linfa/plot_res.py --folder results/ --name TP15_no_disc_error_estimation --iter 10000 --picformat png +#python3 linfa/plot_res.py --folder results/ --name TP15_no_disc_error_estimation --iter 10000 --picformat png # python3 linfa/plot_disc.py --folder results/ --name test_08_lf_w_disc_TP1_uniform_prior --iter 25000 --picformat png --mode histograms --num_points 10 --limfactor 1.0 --saveinterval 1000 --dropouts 10 # python3 linfa/plot_disc.py --folder results/ --name test_19_lf_w_disc_TP15_rep_meas_dropout --iter 10000 --picformat png --mode discr_surface --num_points 10 --limfactor 1.0 --saveinterval 1000 # python3 linfa/plot_disc.py --folder results/ --name test_08_lf_w_disc_TP1_uniform_prior --iter 25000 --picformat png --mode marginal_stats --num_points 10 --limfactor 1.0 --saveinterval 1000 # python3 linfa/plot_disc.py --folder results/ --name TP1_no_disc_gaussian_prior --iter 10000 --picformat png --mode marginal_posterior --num_points 10 --limfactor 1.0 --saveinterval 1000 - +python3 plot_mcmc_and_linfa.py --folder results/ --name TP1_no_disc_uniform_prior --iter 10000 --picformat png