add local analysis test for uncertainty

Author: sevyharris

test/rmgpy/tools/chemDir/solver/sensitivity_1_SPC_18.csv

@@ -0,0 +1,6 @@
+Time (s),dln[C2H6(18)]/dln[k1]: O(0)+H2O2(3)<=>OH(1)+HO2(2),dln[C2H6(18)]/dln[k2]: CH3(14)+PC3H7(15)<=>CH4(16)+CH2CH2CH2(17),dln[C2H6(18)]/dln[k3]: C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19),dln[C2H6(18)]/dln[k4]: C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15),dln[C2H6(18)]/dln[k5]: CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15),dln[C2H6(18)]/dln[k6]: CH3(14)+C2H6(18)<=>CH4(16)+C2H5(12),dln[C2H6(18)]/dln[k7]: HCCO(10)(+M)<=>O(0)+C2H(8)(+M),dln[C2H6(18)]/dln[k8]: C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22),dln[C2H6(18)]/dln[k9]: C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22),dln[C2H6(18)]/dln[k10]: CH3(14)+C2H5(12)<=>CH4(16)+C2H4(11),dln[C2H6(18)]/dG[C2H6(18)],dln[C2H6(18)]/dG[CH4(16)],dln[C2H6(18)]/dG[CH2(5)],dln[C2H6(18)]/dG[CH(4)],dln[C2H6(18)]/dG[C2H3(20)],dln[C2H6(18)]/dG[CH3(14)],dln[C2H6(18)]/dG[C2H4(11)],dln[C2H6(18)]/dG[C2H5(12)]
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+2.74888478006841E-18,-2.54948202221999E-16,-2.09986641193727E-28,-1.68177517988414E-18,9.47815581205925E-43,-4.74240993812363E-29,1.56391028678446E-54,6.05624693198534E-47,-1.10008312605422E-28,3.64619470520282E-56,-3.31477117062407E-66,-9.93392859924497E-17,-9.89301193768683E-48,6.51000897308177E-19,6.51000897344892E-19,1.83574663375722E-29,-4.06622729584823E-53,3.06108016791651E-49,4.22004499424344E-59
+1E-12,-9.27460693505743E-11,-2.0460917755381E-17,-6.1181805764998E-13,3.40947333001306E-26,-4.62135447625443E-18,1.34821506115326E-32,2.66263273777064E-30,-1.07185782928551E-17,3.14317510387031E-34,-7.91446790404999E-39,-3.61380457995889E-11,-2.37250016862895E-31,2.36827788154325E-13,2.36831365928509E-13,1.78888709197875E-18,-8.03285951906872E-37,7.34152354170913E-33,2.43214792225303E-37
+6.14209255203048E-06,-0.0108552687245762,-0.00977366113631576,0.000596136582961715,0.00200221535342703,-0.00846755042343432,0.00013450902403258,0.000108628449374283,-0.00628567915883179,8.65139619102715E-05,-1.74774283067548E-05,-0.00452281325424481,1.26216421675667E-05,-0.00328956503163016,0.00347689688022537,0.00337489330053554,-8.94698376820376E-06,2.89680320928337E-06,6.86042557778715E-07
+0.000511952119227999,-0.259795045759184,-0.127199037182731,0.201488393268386,0.147536983612322,-0.0820373529691973,0.0227660362837704,0.0633512248907625,-0.692292464491573,0.0214205601994878,-0.0568546779535624,-0.119942334769817,0.0416873316334017,-0.146923796990328,0.184660504897779,0.142566434880232,-0.0247332248569192,0.00807885139388224,0.00380827975950797

test/rmgpy/tools/uncertaintyTest.py

@@ -48,8 +48,10 @@ def setup_class(cls):
         chemkin_file = os.path.join(chem_dir, "chem_annotated.inp")
         spc_dict = os.path.join(chem_dir, "species_dictionary.txt")
 
-        cls.uncertainty = Uncertainty(output_directory="chemDir")
+        cls.uncertainty = Uncertainty(output_directory=os.path.abspath(os.path.join(os.path.dirname(__file__), "chemDir")))
         cls.uncertainty.load_model(chemkin_file, spc_dict)
+        for i in range(len(cls.uncertainty.species_list)):
+            cls.uncertainty.species_list[i].index = i  # local analysis depends on species being indexed
 
         # load database properly
         cls.uncertainty.database = RMGDatabase()
@@ -178,6 +180,117 @@ def test_uncertainty_assignment(self):
             rtol=1e-4
         )
 
+    def test_local_analysis(self):
+        """
+        Test to run uncorrelated and then correlated local_analysis and make sure the results are expected
+        """
+        # these are sorted by the name of the corresponding parameter
+        expected_uncorrelated_total_variance = 1.8329056941266446
+        expected_uncorrelated_thermo_variances = np.array([0.09781627, 0.00041446, 9.953e-05, 0.06186124, 0.17092419, 0.04856985, 0.00306632, 0.00391013])
+        expected_uncorrelated_kinetics_variances = np.array([0.00080811, 0.01687337, 0.02022449, 0.15888459, 0.085189, 0.01588366, 0.00012957, 0.01605351, 1.1311145, 0.0010829])
+        expected_correlated_total_variance = 3.006637059881831
+        expected_correlated_thermo_variances = np.array([0.04573167, 0.00014685, 3.263e-05, 0.03236887, 0.07672388, 0.0013764, 0.0001338, 0.00143231, 0.00028968, 0.01098087, 0.00031352, 0.01747643, 0.04856985, 0.09145902])
+        expected_correlated_kinetics_variances = np.array([0.00011471, 0.11981721, 0.15516713, 2.0212174, 0.28939964, 0.00193506, 0.02838292, 0.0161796, 0.0040449, 0.00136045, 0.00168253, 0.00253735, 0.00253735, 0.00136045, 0.01687337, 0.01605351, 0.00080811, 0.00012957])
+        expected_uncorrelated_thermo_labels = {
+            'dln[C2H6(18)]/dG[C2H3(20)]',
+            'dln[C2H6(18)]/dG[C2H4(11)]',
+            'dln[C2H6(18)]/dG[C2H5(12)]',
+            'dln[C2H6(18)]/dG[C2H6(18)]',
+            'dln[C2H6(18)]/dG[CH(4)]',
+            'dln[C2H6(18)]/dG[CH2(5)]',
+            'dln[C2H6(18)]/dG[CH3(14)]',
+            'dln[C2H6(18)]/dG[CH4(16)]',
+        }
+        expected_uncorrelated_kinetics_labels = {
+            'k10: CH3(14)+C2H5(12)<=>CH4(16)+C2H4(11)',
+            'k1: O(0)+H2O2(3)<=>OH(1)+HO2(2)',
+            'k2: CH3(14)+PC3H7(15)<=>CH4(16)+CH2CH2CH2(17)',
+            'k3: C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19)',
+            'k4: C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15)',
+            'k5: CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15)',
+            'k6: CH3(14)+C2H6(18)<=>CH4(16)+C2H5(12)',
+            'k7: HCCO(10)(+M)<=>O(0)+C2H(8)(+M)',
+            'k8: C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22)',
+            'k9: C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22)',
+        }
+        expected_correlated_thermo_labels = {
+            'Estimation C2H3(20)',
+            'Estimation C2H4(11)',
+            'Estimation C2H5(12)',
+            'Estimation C2H6(18)',
+            'Estimation CH(4)',
+            'Estimation CH3(14)',
+            'Group(group) Cds-CdsHH',
+            'Group(group) Cs-CsHHH',
+            'Group(other) R',
+            'Group(radical) CCJ',
+            'Group(radical) CH3',
+            'Group(radical) CJ3',
+            'Library CH2(5)',
+            'Library CH4(16)',
+        }
+        expected_correlated_kinetics_labels = {
+            'Disproportionation Root_Ext-1R!H-R_N-4R->O_N-Sp-5R!H=1R!H_Ext-4CHNS-R_N-6R!H->S_4CHNS->C_N-Sp-6BrBrBrCCCClClClFFFIIINNNOOOPPPSiSiSi#4C_6BrCClFINOPSi->C_N-1R!H-inRing_N-Sp-6C-4C',
+            'Disproportionation Root_Ext-2R!H-R_2R!H->C_4R->C',
+            'Estimation C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15)',
+            'Estimation C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22)',
+            'Estimation C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19)',
+            'Estimation C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22)',
+            'Estimation CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15)',
+            'Estimation CH3(14)+PC3H7(15)<=>CH4(16)+CH2CH2CH2(17)',
+            'H_Abstraction C/H3/Cs;C_methyl',
+            'H_Abstraction C/H3/Cs\\H2\\Cs|O;Cd_Cd\\H2_rad/Cs',
+            'H_Abstraction C/H3/Cs\\H2\\O;C_methyl',
+            'H_Abstraction C/H3/Cs\\H3;C_rad/H2/Cs\\H3',
+            'H_Abstraction C/H3/Cs\\H3;C_rad/H2/Cs\\H\\Cs\\Cs|O',
+            'H_Abstraction C/H3/Cs\\H3;Cd_Cd\\H2_pri_rad',
+            'Library O(0)+H2O2(3)<=>OH(1)+HO2(2)',
+            'PDep HCCO(10)(+M)<=>O(0)+C2H(8)(+M)',
+            'Training Disproportionation CH3(14)+C2H5(12)<=>CH4(16)+C2H4(11)',
+            'Training H_Abstraction CH3(14)+C2H6(18)<=>CH4(16)+C2H5(12)',
+        }
+
+        sensitive_species = [self.uncertainty.species_list[18]]
+
+        # uncorrelated analysis first
+        self.uncertainty.assign_parameter_uncertainties()
+        output = self.uncertainty.local_analysis(sensitive_species=sensitive_species)
+        total_variance, kinetic_uncertainty, thermo_uncertainty = output[sensitive_species[0]]
+        assert np.isclose(total_variance, expected_uncorrelated_total_variance)
+
+        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by name
+        kinetic_variances = [r[2] for r in kinetic_uncertainty]
+        kinetics_names = [r[0] for r in kinetic_uncertainty]
+        sorted_kinetic_variances = [x for _, x in sorted(zip(kinetics_names, kinetic_variances))]
+        assert np.isclose(sorted_kinetic_variances, expected_uncorrelated_kinetics_variances).all()
+        assert set(kinetics_names) == expected_uncorrelated_kinetics_labels
+
+        thermo_variances = [s[2] for s in thermo_uncertainty]
+        thermo_names = [s[0] for s in thermo_uncertainty]
+        sorted_thermo_variances = [x for _, x in sorted(zip(thermo_names, thermo_variances))]
+        assert np.isclose(sorted_thermo_variances, expected_uncorrelated_thermo_variances).all()
+        assert set(thermo_names) == expected_uncorrelated_thermo_labels
+
+        # now repeat for correlated analysis
+        self.uncertainty.assign_parameter_uncertainties(correlated=True)
+        output = self.uncertainty.local_analysis(sensitive_species=sensitive_species, correlated=True)
+        total_variance, kinetic_uncertainty, thermo_uncertainty = output[sensitive_species[0]]
+        assert np.isclose(total_variance, expected_correlated_total_variance)
+
+        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by name
+        kinetic_variances = [r[2] for r in kinetic_uncertainty]
+        kinetics_names = [r[0] for r in kinetic_uncertainty]
+        sorted_kinetic_variances = [x for _, x in sorted(zip(kinetics_names, kinetic_variances))]
+        assert np.isclose(sorted_kinetic_variances, expected_correlated_kinetics_variances).all()
+        assert set(kinetics_names) == expected_correlated_kinetics_labels
+
+        thermo_variances = [s[2] for s in thermo_uncertainty]
+        thermo_names = [s[0] for s in thermo_uncertainty]
+        sorted_thermo_variances = [x for _, x in sorted(zip(thermo_names, thermo_variances))]
+        assert np.isclose(sorted_thermo_variances, expected_correlated_thermo_variances).all()
+        assert set(thermo_names) == expected_correlated_thermo_labels
+
+
     def test_specific_species_uncertainties(self):
         """
         Test uncertainties for a few specific examples