Optimized CompareSpectra such that a model spectrum is only extracted…

… once per grid point when using the compare_model method, using context managers for opening the HDF5 database throughout the database module, replaced the use of spectres with interp1d in ReadModel
tomasstolker · Jul 26, 2023 · b7ab535 · b7ab535
1 parent bc72048
commit b7ab535
Show file tree

Hide file tree

Showing 4 changed files with 856 additions and 883 deletions.
diff --git a/species/analysis/compare_spectra.py b/species/analysis/compare_spectra.py
@@ -553,26 +553,26 @@ def compare_model(
                                 if model_param[5] is not None:
                                     param_dict[model_param[5]] = coord_5_item
 
+                                model_reader = read_model.ReadModel(model)
+
+                                model_box_full = model_reader.get_data(param_dict)
+
                                 for spec_item in self.spec_name:
                                     obj_spec = self.object.get_spectrum()[spec_item][0]
                                     obj_res = self.object.get_spectrum()[spec_item][3]
 
-                                    wavel_range = (
-                                        0.9 * obj_spec[0, 0],
-                                        1.1 * obj_spec[-1, 0],
-                                    )
+                                    # Smooth model spectrum
 
-                                    model_reader = read_model.ReadModel(
-                                        model, wavel_range=wavel_range
+                                    model_flux = read_util.smooth_spectrum(
+                                        model_box_full.wavelength, model_box_full.flux, obj_res
                                     )
 
-                                    model_box = model_reader.get_data(
-                                        param_dict,
-                                        spec_res=obj_res,
-                                        wavel_resample=obj_spec[:, 0],
-                                    )
+                                    # Resample model spectrum
+
+                                    flux_intep = interp1d(model_box_full.wavelength, model_flux)
+                                    model_flux = flux_intep(obj_spec[:, 0])
 
-                                    nan_wavel = np.sum(np.isnan(model_box.flux))
+                                    nan_wavel = np.sum(np.isnan(model_flux))
 
                                     if nan_wavel > 0:
                                         warnings.warn(
@@ -585,7 +585,7 @@ def compare_model(
                                             "calculating the goodness-of-fit statistic."
                                         )
 
-                                    model_spec[spec_item] = model_box.flux
+                                    model_spec[spec_item] = model_flux
 
                                 model_list = []
                                 data_list = []