added AIChE edits, created pairplots

Author: kylajones

aiche_delete_l8r.png

@@ -0,0 +1,174 @@
+import matplotlib.pyplot as plt
+from scipy import stats
+import numpy as np
+
+
+# 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']    = 8
+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'
+
+
+if False:
+    U = np.linspace(-4,4, 200)
+    u = stats.norm.pdf(U, loc = 0, scale = 1)
+    pi = 0.2
+    z = pi * stats.norm.pdf(U, loc = -2, scale = 0.5) + (1 - pi) * stats.norm.pdf(U, loc = 1, scale = 0.5)
+
+    plt.figure()
+    plt.plot(U,u, linewidth = 3)
+    plt.fill_between(U, u * 0, u, alpha = 0.5)
+    plt.savefig("aiche_delete_l8r_blue", dpi = 300)
+    plt.close()
+
+    plt.figure()
+    plt.plot(U,z, linewidth = 3, color = "r")
+    plt.fill_between(U, u * 0, z, alpha = 0.5, color = "r")
+    plt.savefig("aiche_delete_l8r_orange", dpi = 300)
+    plt.close()
+
+if False:
+    X = np.linspace(-5, 5, 200)
+    shift = 6
+    scale = 1
+    x1 = stats.norm.pdf(X, 0 , 1)
+    x2 = 0.75 * stats.norm.pdf(X - shift, 0 - shift, 1.5 * scale)
+    x3 = 0.9 * stats.norm.pdf(X + shift, 0 + shift, 1.25 * scale)
+
+    plt.figure(figsize=(10, 4))
+
+    
+    plt.plot(X - shift, x2, linewidth = 3, color = "r")
+    plt.fill_between(X - shift, x2 * 0, x2, alpha = 0.3, color = "r")
+
+    plt.plot(X + shift, x3, linewidth = 3, color = "b")
+    plt.fill_between(X + shift, x3 * 0, x3, alpha = 0.3, color = "b")
+
+    plt.plot(X, x1, linewidth = 3, color = "g")
+    plt.fill_between(X, x1 * 0, x1, alpha = 0.3, color = "g")
+
+    plt.plot(X, x1*0, linewidth = 3, color = "k")
+    plt.plot(X - shift, x1*0, linewidth = 3, color = "k")
+    plt.plot(X + shift, x1*0, linewidth = 3, color = "k")
+    plt.savefig("kl_divergence", dpi = 300)
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from scipy.stats import norm
+
+    # Observed data
+    np.random.seed(0)
+    x_observed = np.random.normal(loc=2.5, scale=1.0, size=10)  # Data with unknown mean and known variance
+
+    # Known parameters
+    sigma = 1.0  # Known standard deviation of the observations
+    n = len(x_observed)  # Number of observations
+    sample_mean = np.mean(x_observed)
+
+    # Prior parameters
+    mu_0 = 0.0  # Prior mean
+    tau_0 = 1.0  # Prior standard deviation
+
+    # Posterior parameters
+    tau_n_sq = 1 / (1 / tau_0**2 + n / sigma**2)  # Posterior variance
+    mu_n = tau_n_sq * (mu_0 / tau_0**2 + n * sample_mean / sigma**2)  # Posterior mean
+    tau_n = np.sqrt(tau_n_sq)  # Posterior standard deviation
+
+    # Range of mu values
+    mu_values = np.linspace(1.5, 4.5, 200)
+
+    # Calculate unnormalized posterior (prior * likelihood)
+    prior = norm.pdf(mu_values, mu_0, tau_0)
+    likelihood = norm.pdf(mu_values, sample_mean, sigma / np.sqrt(n))
+    unnormalized_posterior = prior * likelihood
+
+    # Scale up unnormalized posterior for visibility
+    unnormalized_posterior_scaled = unnormalized_posterior * 200  # Adjust this factor as needed
+
+    # Calculate normalized posterior
+    normalized_posterior = norm.pdf(mu_values, mu_n, tau_n)
+
+    # Plotting
+    plt.figure(figsize=(5, 5))
+
+    # Unnormalized posterior with shading (scaled up)
+    plt.plot(mu_values, unnormalized_posterior_scaled, label="Unnormalized", linestyle="--", color="blue")
+    plt.fill_between(mu_values, unnormalized_posterior_scaled, color="blue", alpha=0.2)
+
+    # Normalized posterior with shading
+    plt.plot(mu_values, normalized_posterior, label="Normalized", color="red")
+    plt.fill_between(mu_values, normalized_posterior, color="red", alpha=0.2)
+
+    # Annotation to highlight normalization importance
+    plt.text(4.0, 0.8, "Area = 1", color="red", fontsize=16, ha="center", backgroundcolor="white")
+    plt.text(4.0, 0.5, "Area < 1", color="blue", fontsize=16, ha="center", backgroundcolor="white")
+
+    # Labels and legend
+    plt.xlabel(r"$\mu$")
+    plt.ylabel("Density")
+    plt.legend(loc="upper right", fontsize = 14)
+    plt.ylim(0,1.4)
+
+    plt.savefig("normalized_density_fig", dpi = 300)
+
+import os
+import torch
+from linfa.models.discrepancy_models import PhysChem
+import matplotlib.pyplot as plt
+
+
+
+samples = np.loadtxt('results/TP15_no_disc_error_estimation_aiche/TP15_no_disc_error_estimation_aiche_outputs_lf+noise_6000')
+observations = np.loadtxt("observations.csv", skiprows=1, delimiter=',')
+
+
+# Set variable grid
+T = [50.0, 400.0, 450.0]
+P = [1.0, 2.0, 3.0, 4.0, 5.0]
+
+samples = samples.reshape(3,5,5000)
+
+
+# Plot the samples
+for i in range(5000):
+    plt.plot(P, samples[0, :, i], 'm-', linewidth=0.005)
+    plt.plot(P, samples[1, :, i], 'r-', linewidth=0.005)
+    plt.plot(P, samples[2, :, i], color="orange", linestyle='-', linewidth=0.005)
+
+# Plot observations
+plt.plot(observations[:, 1], observations[:, 2], 'ko')
+
+# Add custom legend entries by plotting representative lines
+plt.plot([], [], 'm-', label="350 K")       # Magenta line
+plt.plot([], [], 'r-', label="400 K")       # Red line
+plt.plot([], [], color="orange", linestyle='-', label="450 K")  # Orange line
+plt.plot([], [], 'ko', label="Observations")  # Black circles for observations
+
+# Add legend, labels, and limits
+plt.legend()
+plt.xlim(1, 5)
+plt.xlabel("Pressure, $P$ [Pa]")
+plt.ylabel("Coverage")
+plt.savefig("function", dpi=300)
+
+# for i in range(5000):
+#     plt.plot(P, samples[0,:,i],'m-', linewidth = 0.005)
+#     plt.plot(P, samples[1,:,i],'r-', linewidth = 0.005)
+#     plt.plot(P, samples[2,:,i],color = "orange", linestyle = '-', linewidth = 0.005)
+# plt.plot(observations[:,1], observations[:,2], 'ko')
+# plt.xlim(1,5)
+# plt.xlabel("Pressure, $P$ [Pa]")
+# plt.ylabel("Coverage")
+# plt.savefig("function", dpi = 300)
+# exit()

aiche_delete_l8r_blue.png

@@ -0,0 +1,61 @@
+# import numpy as np
+# import os
+
+# # Define file path
+# filename = "TP15_no_disc_error_estimation_2"
+# path = os.path.join("results", filename)
+# iters = 10000
+
+# # Load data
+# mcmc = np.loadtxt(os.path.join(path, 'mcmc'))
+# linfa = np.loadtxt(os.path.join(path, f"{filename}_params_{iters}"))
+
+# # Ground truth parameters
+# gt_params = np.array([1.0E3, -21.0E3, 0.05])
+
+# # Calculate RMSE for mcmc
+# diff_mcmc = (mcmc - gt_params) ** 2
+# rmse_mcmc = np.sqrt(np.mean(diff_mcmc))
+
+# # Calculate RMSE for linfa
+# diff_linfa = (linfa - gt_params) ** 2
+# rmse_linfa = np.sqrt(np.mean(diff_linfa))
+
+# print("RMSE for MCMC:", rmse_mcmc)
+# print("RMSE for Linfa:", rmse_linfa)
+
+import numpy as np
+import os
+
+# Define file path
+filename = "TP15_no_disc_error_estimation_aiche"
+path = os.path.join("results", filename)
+iters = 6000
+
+# Load data
+mcmc = np.loadtxt(os.path.join(path, 'mcmc'))
+linfa = np.loadtxt(os.path.join(path, f"{filename}_params_{iters}"))
+
+# Ground truth parameters
+gt_params = np.array([1.0E3, -21.0E3, 0.05])
+
+# Calculate RMSE for mcmc
+diff_mcmc = (mcmc - gt_params) ** 2
+rmse_mcmc = np.sqrt(np.mean(diff_mcmc))
+
+# Calculate RMSE for linfa
+diff_linfa = (linfa - gt_params) ** 2
+rmse_linfa = np.sqrt(np.mean(diff_linfa))
+
+# Calculate the average standard deviation of β samples for mcmc
+std_dev_mcmc = np.std(mcmc, axis=0)  # Standard deviation for each parameter
+avg_beta_sd_mcmc = np.mean(std_dev_mcmc)  # Average standard deviation across parameters
+
+# Calculate the average standard deviation of β samples for linfa
+std_dev_linfa = np.std(linfa, axis=0)  # Standard deviation for each parameter
+avg_beta_sd_linfa = np.mean(std_dev_linfa)  # Average standard deviation across parameters
+
+print("RMSE for MCMC:", rmse_mcmc)
+print("RMSE for Linfa:", rmse_linfa)
+print("Average standard deviation of β samples for MCMC:", avg_beta_sd_mcmc)
+print("Average standard deviation of β samples for Linfa:", avg_beta_sd_linfa)

aiche_delete_l8r_orange.png

@@ -119,10 +119,11 @@ def solve_true(self, cal_inputs):
 
         # compute equilibrium constant of site one
         k1Const = 1.0/p01Const * torch.exp(-e1Const / self.RConst / T)
+        print(k1Const)
 
         # compute equilibrium constant of site two
         k2Const = 1.0/p02Const * torch.exp(-e2Const / self.RConst / T)
-
+        print(k2Const)
         # compute surface coverage fraction for two adsorption sites with different equilibrium conditions
         cov_frac = lambda1Const * (k1Const*P/(1 + k1Const*P)) + lambda2Const * (k2Const*P/(1 + k2Const*P))
 

aiche_plots.py

@@ -26,15 +26,18 @@ def plot_log(log_file,out_dir,fig_format='png', use_dark_mode=False):
   log_data = np.loadtxt(log_file)
   # loss profile
   plt.figure()
-  plt.plot(log_data[:, 1],log_data[:, 2],'b-', label = 'ELBO')
-  plt.plot(log_data[:, 1],log_data[:, 3],'g-', label = 'Sum Log Jacobian')
-  plt.plot(log_data[:, 1],log_data[:, 4],'m-', label = 'Likelihood')
-  plt.plot(log_data[:, 1],log_data[:, 5],'k-', label = 'Prior')
+
+  plt.plot(log_data[:, 1], log_data[:, 2],'b-', label = 'ELBO')
+  # plt.plot(log_data[:, 1],log_data[:, 3],'g-', label = 'Sum Log Jacobian')
+  # plt.plot(log_data[:, 1],log_data[:, 4],'m-', label = 'Likelihood')
+  # plt.plot(log_data[:, 1],log_data[:, 5],'k-', label = 'Prior')
   plt.xlabel('Iterations')
   plt.ylabel('Log Loss')
   plt.legend()
   plt.savefig(os.path.join(out_dir,'log_plot.'+fig_format))
   plt.close()
+  print("success")
+  exit()
 
 def plot_params(param_data,LL_data,idx1,idx2,out_dir,out_info,fig_format='png', use_dark_mode=False):  
 

error_calcs.py

@@ -13,10 +13,10 @@
 def run_test():
 
     exp = experiment()
-    exp.name = "TP15_no_disc_error_estimation"
+    exp.name = "TP15_no_disc_error_estimation_2"
     exp.flow_type           = 'realnvp'     # str: Type of flow (default 'realnvp') # TODO: generalize to work for TP1
-    exp.n_blocks            = 50            # int: Number of hidden layers  
-    exp.hidden_size         = 10            # int: Hidden layer size for MADE in each layer (default 100)
+    exp.n_blocks            = 1            # int: Number of hidden layers  
+    exp.hidden_size         = 10           # int: Hidden layer size for MADE in each layer (default 100)
     exp.n_hidden            = 1             # int: Number of hidden layers in each MADE
     exp.activation_fn       = 'relu'        # str: Activation function used (default 'relu')
     exp.input_order         = 'sequential'  # str: Input oder for create_mask (default 'sequential')
@@ -25,10 +25,10 @@ def run_test():
 
     # p0,e,sigma_e (measurement noise also estimated)
     exp.input_size          = 3             # int: Dimensionalty of input (default 2)
-    exp.batch_size          = 300           # int: Number of samples generated (default 100)
+    exp.batch_size          = 100           # int: Number of samples generated (default 100)
     exp.true_data_num       = 1             # double: Number of true model evaluted (default 2)
-    exp.n_iter              = 4000          # int: Number of iterations (default 25001)
-    exp.lr                  = 0.0005        # float: Learning rate (default 0.003)
+    exp.n_iter              = 10000          # int: Number of iterations (default 25001)
+    exp.lr                  = 0.001        # float: Learning rate (default 0.003)
     exp.lr_decay            = 0.9999        # float:  Learning rate decay (default 0.9999)
     exp.log_interval        = 10            # int: How often to show loss stat (default 10)
 
@@ -56,9 +56,10 @@ def run_test():
     exp.device = torch.device('cuda:0' if torch.cuda.is_available() and not exp.no_cuda else 'cpu')
 
     # Define transformation
-    trsf_info = [['tanh', -40.0, 30.0, 500.0, 5000.0],
+    trsf_info = [['tanh', -50.0, 100.0, 500.0, 1500.0],
                  ['tanh', -30.0, 30.0, -30000.0, -15000.0],
-                 ['tanh', -20.0, 20.0, 0.00001, 0.2]]
+                 ['tanh', -10.0, 20.0, 0.00001, 0.2]]
+
     trsf = Transformation(trsf_info)
 
     # Apply the transformation
@@ -152,14 +153,12 @@ def log_prior(calib_inputs, transform):
         # Compute the calibration inputs in the physical domain
         phys_inputs = transform.forward(calib_inputs)
 
-        # Gaussian prior on physically meaningful inputs
-        mean = np.array([1.0E3, -21.0E3])
-        cov = np.array([[100.0**2, 0.0],
-                        [0.0,   500.0**2]])
-        gauss_prior_res = torch.from_numpy(stats.multivariate_normal.logpdf(phys_inputs[:,:2].detach().numpy(), mean = mean, cov = cov))
-        
+        # Gaussian prior on physically meaningful inputs        
+        mean = torch.tensor([1.0E3, -21.0E3])
+        cov = torch.diag(torch.tensor([100.0**2, 500.0**2]))
+        gauss_prior_res = torch.distributions.MultivariateNormal(mean, cov).log_prob(phys_inputs[:,:2])
+
         # Beta prior on noise
-        #torch.distributions.gamma.Gamma(torch.tensor([2.0]), torch.tensor([2.0]))
         std_dev_ratio = torch.distributions.beta.Beta(torch.tensor([1.0]), torch.tensor([19.0]))
 
         if False:
@@ -196,6 +195,6 @@ def generate_data(use_true_model = False, num_observations=50):
 # Main code
 if __name__ == "__main__":
 
-    #generate_data(use_true_model = False, num_observations = 1)
+    # generate_data(use_true_model = False, num_observations = 1)
 
     run_test()