generalized plotting to work with p parameters

linfa/tests/test_no_disc_TP15_error_estimation.py

@@ -131,11 +131,6 @@ def log_density(calib_inputs, model, surrogate, transform):
             # - 1 / (2 * sigma^2) sum_{i = 1} ^ N (eta_i + disc_i - y_i)^2 
             l3 = -0.5 / (std_dev ** 2) * torch.sum((modelOut + discrepancy.t() - Data[:,loopA].unsqueeze(0))**2, dim = 1)
 
-            # Compute negative ll (num_batch x 1)
-            # print('l1', l1)
-            # print('l2', l2)
-            # print('l3', l3)
-            # exit()
             negLL = -(l1 + l2 + l3) # sum contributions
             res = -negLL.reshape(calib_inputs.size(0), 1) # reshape
 
@@ -180,9 +175,9 @@ def log_prior(calib_inputs, transform):
         # Add gaussian log prior for first two parameters
         gauss_prior_res = l1 + l2 + l3
 
-        if False:
+        if True:
              # Add beta prior for third parameter
-            sigma_prior = torch.distributions.beta.Beta(torch.tensor([1.0]), torch.tensor([0.5]))
+            sigma_prior = torch.distributions.beta.Beta(torch.tensor([3/16]), torch.tensor([57/16]))
             prior_res = sigma_prior.log_prob(phys_inputs[:,2])
         else:
             # Add uniform prior

plot_mcmc_and_linfa.py

@@ -22,13 +22,30 @@
 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):  
+def plot_marginals(param_data, idx1, fig_format='png'):
+
+  param_data = np.loadtxt(param_data)  # [5000, 3]
+  mcmc_1_data = np.loadtxt(os.path.join(out_dir,'posterior_samples_'+str(idx1)))
+  gt_params = [1000, -21.0E3, 0.05]
+
+  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(gt_params[idx1], 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)+'.'+fig_format))
+  plt.close()
+
+
+def plot_params(param_data, LL_data, idx1, idx2, 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)
+  mcmc_1_data = np.loadtxt(os.path.join(out_dir,'posterior_samples_'+str(idx1)))
+  mcmc_2_data = np.loadtxt(os.path.join(out_dir,'posterior_samples_'+str(idx2)))
+  gt_params = [1000, -21.0E3, 0.05]
 
   # Combine MCMC samples
   samples = np.vstack([mcmc_1_data, mcmc_2_data])  # Transpose to get shape (n, d)
@@ -47,31 +64,13 @@ def plot_params(param_data, LL_data, idx1, idx2, param_file_mcmc_1, param_file_m
   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.plot(gt_params[idx1], gt_params[idx2], '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
 # =========
@@ -142,16 +141,15 @@ def plot_params(param_data, LL_data, idx1, idx2, param_file_mcmc_1, param_file_m
   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
+      plot_marginals(param_file, loopA)
       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)
+        plot_params(param_file,LL_file,loopA,loopB,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)

run_plot_res.sh

@@ -3,5 +3,5 @@
 # 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
+python3 plot_mcmc_and_linfa.py --folder results/ --name TP15_no_disc_error_estimation --iter 10000 --picformat png