Merge assign_parameter_uncertainties and intermediate_uncertainties

There's a lot of repeated code between the old UQ's assign_parameter_uncertainties and the new UQ's assign_intermediate_uncertainties, so this combines them into a single function in order to reuse a lot of that repeated code and make it easier to keep them consistent.

Author: sevyharris

rmgpy/tools/uncertainty.py

@@ -916,43 +916,56 @@ def assign_parameter_uncertainties(self, g_param_engine=None, k_param_engine=Non
         self.thermo_input_uncertainties = []
         self.kinetic_input_uncertainties = []
 
+        self.dG_dqs = []  # a list of dictionaries to store the intermediate derivatives dG_i/dq for each parameter q that contributes to the uncertainty of species i's G
+        self.dlnk_dqs = []  # a list of dictionaries to store the intermediate derivatives dlnk_i/dq for each parameter q that contributes to the uncertainty of reaction i's k
+
         for species in self.species_list:
             if not correlated:
                 dG = g_param_engine.get_uncertainty_value(self.species_sources_dict[species])
                 self.thermo_input_uncertainties.append(dG)
             else:
                 source = self.species_sources_dict[species]
                 dG = {}
+                dG_dq = {}
                 if 'Library' in source:
                     pdG = g_param_engine.get_partial_uncertainty_value(source, 'Library', corr_param=source['Library'])
                     label = source['Library'][2]  # the label we added to the end of the library source tuple in extract_sources_from_model
                     dG[label] = pdG
+                    dG_dq[label] = 1  # derivative of G with respect to the library parameter is 1, since it's a direct contribution
                 if 'Surface_Library' in source:
                     pdG = g_param_engine.get_partial_uncertainty_value(source, 'Surface_Library', corr_param=source['Surface_Library'])
                     label = source['Surface_Library'][2]
                     dG[label] = pdG
+                    dG_dq[label] = 1  # derivative of G with respect to the surface library parameter is 1, since it's a direct contribution
                 if 'QM' in source:
                     pdG = g_param_engine.get_partial_uncertainty_value(source, 'QM', corr_param=source['QM'])
                     label = source['QM'][2]
                     dG[label] = pdG
+                    dG_dq[label] = 1  # derivative of G with respect to the QM parameter is 1, since it's a direct contribution
                 if 'ADS' in source:
                     for adsGroupType, groupList in source['ADS'].items():
                         for group, weight in groupList:
                             pdG = g_param_engine.get_partial_uncertainty_value(source, 'ADS', group, adsGroupType)
                             label = 'AdsorptionCorrection({}) {}'.format(adsGroupType, group.label)
                             dG[label] = pdG
+                            if weight != 1:
+                                raise ValueError('Weight for adsorption group contribution to thermo should be 1, but got weight={weight} for {adsGroupType} in species {species}'.format(weight=weight, adsGroupType=adsGroupType, species=species))
+                            dG_dq[label] = weight  # This should be 1 because there's only one group contribution per adsorption correction
                 if 'GAV' in source:
                     for groupType, groupList in source['GAV'].items():
                         for group, weight in groupList:
                             pdG = g_param_engine.get_partial_uncertainty_value(source, 'GAV', group, groupType)
                             label = 'Group({}) {}'.format(groupType, group.label)
                             dG[label] = pdG
+                            dG_dq[label] = weight  # the derivative of G with respect to the group contribution is the weight of that group contribution to the overall G
                     # We also know if there is group additivity used, there will be uncorrelated estimation error
                     est_pdG = g_param_engine.get_partial_uncertainty_value(source, 'Estimation')
                     if est_pdG:
                         label = 'Estimation {}'.format(species.to_chemkin())
                         dG[label] = est_pdG
+                        dG_dq[label] = 1  # the derivative of G with respect to the estimation error is 1, since we add the term directly
                 self.thermo_input_uncertainties.append(dG)
+                self.dG_dqs.append(dG_dq)
 
         for reaction in self.reaction_list:
             if not correlated:
@@ -961,6 +974,7 @@ def assign_parameter_uncertainties(self, g_param_engine=None, k_param_engine=Non
             else:
                 source = self.reaction_sources_dict[reaction]
                 dlnk = {}
+                dlnk_dq = {}
                 if 'Rate Rules' in source:
                     family = source['Rate Rules'][0]
                     source_dict = source['Rate Rules'][1]
@@ -975,166 +989,55 @@ def assign_parameter_uncertainties(self, g_param_engine=None, k_param_engine=Non
                                                                              corr_family=family)
                         label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)
                         dlnk[label] = dplnk
-
+                        dlnk_dq[label] = weight  # the derivative of ln(k) with respect to the rate rule contribution
                     for ruleEntry, trainingEntry, weight in training:
                         dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Rate Rules', corr_param=ruleEntry,
                                                                              corr_family=family)
                         label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)
                         dlnk[label] = dplnk
+                        dlnk_dq[label] = weight  # the derivative of ln(k) with respect to the training rule contribution
+
+                    N = len(source_dict['rules']) + len(source_dict['training'])
 
                     # There is also estimation error if rate rules are used (nonexact and family contribute to this)
                     nonexact_dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Estimation Nonexact', corr_family=family)
                     if nonexact_dplnk:
                         label = 'Estimation Nonexact {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
                         dlnk[label] = nonexact_dplnk
+                        dlnk_dq[label] = np.log10(N + 1)  # the derivative of ln(k) with respect to the nonexact estimation error
 
                     family_dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Estimation Family', corr_family=family)
                     if family_dplnk:
                         label = 'Estimation Family {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
                         dlnk[label] = family_dplnk
-
+                        dlnk_dq[label] = 1  # the derivative of ln(k) with respect to the family estimation error
                 elif 'PDep' in source:
                     dplnk = k_param_engine.get_partial_uncertainty_value(source, 'PDep', source['PDep'])
                     label = 'PDep {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
                     dlnk[label] = dplnk
+                    dlnk_dq[label] = 1
 
                 elif 'Library' in source:
                     dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Library', source['Library'])
                     label = source['Library'][2]
                     dlnk[label] = dplnk
+                    dlnk_dq[label] = 1
 
                 elif 'Surface_Library' in source:
                     dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Surface_Library', source['Surface_Library'])
                     label = source['Surface_Library'][2]
                     dlnk[label] = dplnk
+                    dlnk_dq[label] = 1
 
                 elif 'Training' in source:
                     dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Training', source['Training'])
                     family = source['Training'][0]
                     label = 'Training {} {}'.format(family, reaction.to_chemkin(self.species_list, kinetics=False))
                     dlnk[label] = dplnk
+                    dlnk_dq[label] = 1
 
                 self.kinetic_input_uncertainties.append(dlnk)
-
-    def assign_intermediate_uncertainties(self, g_param_engine=None, k_param_engine=None, correlated=False):
-        """
-        Assign uncertainties to the intermediate parameters based on the sources of the species thermo and reaction kinetics.
-
-        This fills out the class variables thermo_intermediate_uncertainties and kinetic_intermediate_uncertainties
-        these are each list of dictionaries. For every species or reaction, it lists all the intermediate sources contributing to that parameter's uncertainty.
-
-        So for example, thermo_intermediate_uncertainties might look something like this:
-        
-        thermo_intermediate_uncertainties = [
-            {'Group(group) Cds-CdsHH': 2.0, 'Group(radical) CCJ': 1.0, 'Estimation CH(4)': 1.0},
-            {'Library CH2(5)': 1.0},
-        ]
-        The keys of the dictionaries are the label names for the intermediate parameters.
-        and the values are partial derivatives dG_i/dq_w, how the species i Gibbs uncertainty changes with the intermediate parameter w.
-        
-        This function is the new formulation's equivalent to assign_parameter_uncertainties and similarly handles both correlated and uncorrelated cases.
-        But instead of assuming all underlying parameters are independent, here we can allow for dependence as long as we have the covariance
-        """
-        if g_param_engine is None:
-            g_param_engine = ThermoParameterUncertainty(other_covariances=self.thermo_covariances_dict)
-        if k_param_engine is None:
-            k_param_engine = KineticParameterUncertainty()
-
-        self.thermo_intermediate_uncertainties = []  # store the intermediate dG_i/dq for each parameter q that contributes to the uncertainty of G_i, for use in correlated uncertainty analysis
-        self.kinetic_intermediate_uncertainties = []
-
-        for species in self.species_list:
-            if not correlated:
-                dG = g_param_engine.get_uncertainty_value(self.species_sources_dict[species])
-                self.thermo_intermediate_uncertainties.append(dG)  # in the uncorrelated case, the intermediate is just the uncertainty value itself, since there is only one parameter that contributes to the uncertainty
-            else:
-                source = self.species_sources_dict[species]
-                dGdq = {}
-                if 'Library' in source:
-                    label = source['Library'][2]
-                    dGdq[label] = 1  # dG/dG_lib = 1, because the parameter is never scaled by anything other than 1 when it is used
-                if 'Surface_Library' in source:
-                    label = source['Surface_Library'][2]
-                    dGdq[label] = 1  # dG/dG_surf = 1, because the parameter is never scaled by anything other than 1 when it is used
-                if 'QM' in source:
-                    label = source['QM'][2]
-                    dGdq[label] = 1
-                if 'ADS' in source:
-                    for adsGroupType, groupList in source['ADS'].items():
-                        for group, weight in groupList:
-                            label = 'AdsorptionCorrection({}) {}'.format(adsGroupType, group.label)
-                            if weight != 1:
-                                raise ValueError('Weight for adsorption group contribution to thermo should be 1, but got weight={weight} for {adsGroupType} in species {species}'.format(weight=weight, adsGroupType=adsGroupType, species=species))
-                            dGdq[label] = weight  # This should be 1
-                if 'GAV' in source:
-                    for groupType, groupList in source['GAV'].items():
-                        for group, weight in groupList:
-                            label = 'Group({}) {}'.format(groupType, group.label)
-                            dGdq[label] = weight  # dG/dG_group = weight, because the group contribution is scaled by the weight when it is used in the thermo estimation
-                    # We also know if there is group additivity used, there will be uncorrelated estimation error
-                    label = 'Estimation {}'.format(species.to_chemkin())
-                    dGdq[label] = 1  # dG/dG_est = 1, because the estimation error is added on top of the group additivity value, so it is never scaled by anything other than 1 when it is used
-
-                self.thermo_intermediate_uncertainties.append(dGdq)
-
-        for reaction in self.reaction_list:
-            if not correlated:
-                dlnk = k_param_engine.get_uncertainty_value(self.reaction_sources_dict[reaction])
-                self.kinetic_intermediate_uncertainties.append(dlnk)  # in the uncorrelated case, the intermediate is just the uncertainty value itself, since there is only one parameter that contributes to the uncertainty
-            else:
-                source = self.reaction_sources_dict[reaction]
-                dlnkdq = {}
-                if 'Rate Rules' in source:
-                    family = source['Rate Rules'][0]
-                    source_dict = source['Rate Rules'][1]
-                    rules = source_dict['rules']
-                    training = source_dict['training']
-                    exact = source_dict['exact']
-                    surface_prefix = ''
-                    if reaction.is_surface_reaction():
-                        surface_prefix = 'Surface '
-                    for ruleEntry, weight in rules:
-                        label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)
-                        dlnkdq[label] = weight  # dlnk/dlnk_rule = weight, because the rate rule is scaled by the weight when it is used in the kinetics estimation
-
-                    for ruleEntry, trainingEntry, weight in training:
-                        # TODO - test that training reactions in a tree are correlated with the exact match kind of training reaction
-                        # for now, we follow the old convention of treating these as rate rules
-                        label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)  # ruleEntry should probably be the reaction equation itself
-                        dlnkdq[label] = weight  # dlnk/dlnk_training = weight, because the training entry is scaled by the weight when it is used in the kinetics estimation
-
-                    # There is also estimation error if rate rules are used
-                    # Record dlnk/dlnk_family, the derivative with respect to the family estimation uncertainty
-                    label = 'Estimation Family {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                    dlnkdq[label] = 1  # dlnk/dlnk_family = 1, because the family estimation uncertainty is added on top of the rate rule values, so it is never scaled by anything other than 1 when it is used
-
-                    # Record the non-exact estimation error if not an exact match for a rate rule
-                    if not exact:
-                        N = len(source_dict['rules']) + len(source_dict['training'])
-                        label = 'Estimation Nonexact {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                        dlnkdq[label] = np.log10(N + 1) 
-
-                elif 'PDep' in source:
-                    label = 'PDep {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                    dlnkdq[label] = 1.0  # dlnk/dlnk_PDep = 1, because the PDep kinetics is never scaled by anything other than 1 when it is used
-
-                elif 'Library' in source:
-                    label = 'Library {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                    dlnkdq[label] = 1.0  # dlnk/dlnk_lib = 1, because the library kinetics is never scaled by anything other than 1 when it is used
-
-                elif 'Surface_Library' in source:
-                    label = 'Surface_Library {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                    dlnkdq[label] = 1.0  # dlnk/dlnk_surf_lib = 1, because the surface library kinetics is never scaled by anything other than 1 when it is used
-
-                elif 'Training' in source:
-                    family = source['Training'][0]
-                    surface_prefix = ''
-                    if reaction.is_surface_reaction():
-                        surface_prefix = 'Surface '
-                    label = '{}Training {} {}'.format(surface_prefix, family, reaction.to_chemkin(self.species_list, kinetics=False))
-                    dlnkdq[label] = 1.0
-
-                self.kinetic_intermediate_uncertainties.append(dlnkdq)
+                self.dlnk_dqs.append(dlnk_dq)
 
     def sensitivity_analysis(self, initial_mole_fractions, sensitive_species, T, P, termination_time,
                              sensitivity_threshold=1e-3, number=10, fileformat='.png', initial_surface_coverages=None,