added standard deviation as additional parameter to estimate

linfa/models/discrepancy_models.py

@@ -2,9 +2,60 @@
 import itertools
 import numpy as np
 
-class PhysChem(object):
+class PhysChem_general(object):
 
     def __init__(self, var_inputs):
+        pass
+
+    def solve_t(self, cal_inputs):
+        pass
+
+    def solve_true(self, cal_inputs):
+        pass
+
+    def genDataFile(self, dataFileNamePrefix='observations', use_true_model=True, store=True, num_observations=10):
+
+        # solve model
+        if(use_true_model):
+            # TODO: update call to solve_true to include additional parameters
+            def_out = self.solve_true(self.defParams)
+        else:
+            # TODO: update call to solve_t to include additional parameters
+            def_out = self.solve_t(self.defParams)
+
+        # get standard deviation
+        stdDev = self.stdRatio * torch.mean(torch.abs(def_out))
+
+        # add noise to coverage data
+        coverage = def_out.repeat(1, num_observations)
+
+        for loopA in range(num_observations):
+            noise = torch.randn((len(coverage),1)) * stdDev
+            coverage[:,loopA] = coverage[:,loopA] + noise.flatten()
+
+        if store:
+
+            # specify file name
+            file_name = f'{dataFileNamePrefix}.csv'
+
+            with open(file_name, 'w') as file:
+
+                # write file
+                file.write('temperature, pressure, coverage \n')
+
+                # specify row names
+                for loopA in range(self.var_in.size(0)):
+                    for loopB in range(self.var_in.size(1)):
+                        file.write(f'{self.var_in[loopA,loopB]},')                        
+                    for loopB in range(coverage.size(1)-1):
+                        file.write(f'{coverage[loopA,loopB].item()},')
+                    file.write(f'{coverage[loopA,coverage.size(1)-1].item()}\n')
+
+class PhysChem(PhysChem_general):
+
+    def __init__(self, var_inputs):
+
+        super().__init__(var_inputs)
 
         # variable inputs
         self.var_in = torch.tensor(list(itertools.product(*var_inputs)))
@@ -77,45 +128,84 @@ def solve_true(self, cal_inputs):
 
         # Return
         return cov_frac
-
 
-    def genDataFile(self, dataFileNamePrefix='observations', use_true_model=True, store=True, num_observations=10):
+class PhysChem_error(PhysChem_general):
 
-        # solve model
-        if(use_true_model):
-            # TODO: update call to solve_true to include additional parameters
-            def_out = self.solve_true(self.defParams)
-        else:
-            # TODO: update call to solve_t to include additional parameters
-            def_out = self.solve_t(self.defParams)
+    def __init__(self, var_inputs):
+
+        super().__init__(var_inputs)
 
-        # get standard deviation
-        stdDev = self.stdRatio * torch.mean(torch.abs(def_out))
+        # variable inputs
+        self.var_in = torch.tensor(list(itertools.product(*var_inputs)))
 
-        # add noise to coverage data
-        coverage = def_out.repeat(1, num_observations)
+        # calibration inputs
+        self.defParams = torch.tensor([[1e3, -21e3, 0.05]]) # standard presssure (MPa) and energy of ads. (J/mol)
+        self.limits = [[5e2, 1.5e3],[-30e3, -15e3],[0.01,0.3]] # range on defParams
+
+        ## model constants
+        self.RConst = torch.tensor(8.314) # universal gas constant (J/ mol/ K)
+        self.data = None # dataset of model
+        self.stdRatio = 0.05 # standard deviation ratio
+        self.defOut = self.solve_t(self.defParams)
+
+    def solve_t(self, cal_inputs):
 
-        for loopA in range(num_observations):
-            noise = torch.randn((len(coverage),1)) * stdDev
-            coverage[:,loopA] = coverage[:,loopA] + noise.flatten()
+        # TODO: update cal_inputs to pass other variables for inference
+
+        num_batch = len(cal_inputs)
+        num_vars = len(self.var_in)
 
-        if store:
+        # unpack variable inputs
+        T, P = torch.chunk(self.var_in, chunks = 2, dim = 1)
+        T = T.repeat(1, num_batch)
+        P = P.repeat(1, num_batch)
 
-            # specify file name
-            file_name = f'{dataFileNamePrefix}.csv'
+        # unpack calibration inputs
+        p0Const, eConst, sigma_e = torch.chunk(cal_inputs, chunks = 3, dim = 1) # split cal_inputs into two chunks along the second dimension
+        p0Const = p0Const.repeat(1, num_vars).t()
+        eConst  = eConst.repeat(1, num_vars).t()
+
+        # compute equilibrium constant
+        kConst = 1/p0Const* torch.exp(-eConst / self.RConst / T)
+
+        # compute surface coverage fraction
+        cov_frac = kConst*P/(1 + kConst*P)
 
-            with open(file_name, 'w') as file:
+        # Return coverages        
+        return cov_frac
 
-                # write file
-                file.write('temperature, pressure, coverage \n')
+    # TODO: update cal_inputs to include variable number of parameters not called by the model
+    def solve_true(self, cal_inputs):
 
-                # specify row names
-                for loopA in range(self.var_in.size(0)):
-                    for loopB in range(self.var_in.size(1)):
-                        file.write(f'{self.var_in[loopA,loopB]},')                        
-                    for loopB in range(coverage.size(1)-1):
-                        file.write(f'{coverage[loopA,loopB].item()},')
-                    file.write(f'{coverage[loopA,coverage.size(1)-1].item()}\n')
+        num_batch = len(cal_inputs)
+        num_vars = len(self.var_in)
+
+        # unpack variable input
+        T, P = torch.chunk(self.var_in, chunks = 2, dim = 1)
+        T = T.repeat(1, num_batch)
+        P = P.repeat(1, num_batch)
+
+        # unpack ground truth calibration input (adsorption site 1)
+        p01Const, e1Const, sigma_e = torch.chunk(cal_inputs, chunks = 3, dim = 1)
+        p01Const = p01Const.repeat(1, num_vars).t()
+        e1Const  = e1Const.repeat(1, num_vars).t()
+
+        # specify adsorption site two parameters
+        p02Const, e2Const, lambda1Const, lambda2Const = torch.tensor(5000.0), torch.tensor(-22000.0), torch.tensor(1.0), torch.tensor(0.5)
+        p02Const = p02Const.repeat(1, num_vars).t()
+        e2Const  = e2Const.repeat(1, num_vars).t()
+
+        # compute equilibrium constant of site one
+        k1Const = 1.0/p01Const * torch.exp(-e1Const / self.RConst / T)
+
+        # compute equilibrium constant of site two
+        k2Const = 1.0/p02Const * torch.exp(-e2Const / self.RConst / T)
+
+        # 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))
+
+        # Return
+        return cov_frac
 
 # TEST SURROGATE
 if __name__ == '__main__':

linfa/run_experiment.py

@@ -327,12 +327,20 @@ def train(self, nf, optimizer, iteration, log, sampling=True, t=1):
 
         # Free energy bound
         if(self.model_logprior is None):
+            # test1 = (- torch.sum(sum_log_abs_det_jacobians, 1)).mean()
+            # test2 = (self.model_logdensity(xk)).mean()
+            # print(test1,test2)
             loss = (- torch.sum(sum_log_abs_det_jacobians, 1) - t * self.model_logdensity(xk)).mean()
         else:
+            # test1 = (- torch.sum(sum_log_abs_det_jacobians, 1)).mean()
+            # test2 = (self.model_logdensity(xk)).mean()
+            # test3 = (self.model_logprior(xk)).mean()
+            # print(test1,test2,test3)
             loss = (- torch.sum(sum_log_abs_det_jacobians, 1) - t * (self.model_logdensity(xk) + self.model_logprior(xk))).mean()
         optimizer.zero_grad()
         loss.backward()
         optimizer.step()
+
         if iteration % self.log_interval == 0:
             print('VI NF (t=%5.3f): it: %7d | loss: %8.3e' % (t,iteration, loss.item()))
             log.append([t, iteration, loss.item()])