From 87ebe5601442a90454113966f6edebcd128d9a05 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Mon, 9 Oct 2023 11:02:04 +0200 Subject: [PATCH] Support for only including photometric fluxes when comparing a grid of model spectra with CompareSpectra --- species/analysis/compare_spectra.py | 34 +++++++++++++++++++++++++---- species/data/database.py | 29 +++++++++++++++++------- species/data/model_data.json | 4 ++-- species/plot/plot_comparison.py | 26 ++++++++++++++++++---- species/plot/plot_spectrum.py | 6 ++++- species/util/plot_util.py | 1 + 6 files changed, 81 insertions(+), 19 deletions(-) diff --git a/species/analysis/compare_spectra.py b/species/analysis/compare_spectra.py index 1abbbd72..2f890406 100644 --- a/species/analysis/compare_spectra.py +++ b/species/analysis/compare_spectra.py @@ -444,7 +444,11 @@ def compare_model( else: w_i[spec_item] = np.ones(obj_wavel.shape[0]) + # Open the HDF5 databse + species_db = database.Database() + if isinstance(inc_phot, bool): + # The argument of inc_phot is a boolean if inc_phot: # Select all filters if inc_phot=True species_db = database.Database() @@ -454,6 +458,10 @@ def compare_model( else: inc_phot = [] + else: + # The argument of inc_phot is a list with filters + object_box = species_db.get_object(self.object_name, inc_phot=inc_phot) + if scale_spec is None: scale_spec = [] @@ -592,21 +600,26 @@ def compare_model( model_spec[spec_item] = model_flux + g_fit = 0.0 + model_list = [] data_list = [] weights_list = [] + for spec_item in self.spec_name: if spec_item not in scale_spec: model_list.append(model_spec[spec_item]) data_list.append(spec_data[spec_item][0]) weights_list.append(w_i[spec_item]) + model_phot = {} + for phot_item in inc_phot: syn_phot = photometry.SyntheticPhotometry(phot_item) - model_phot = syn_phot.spectrum_to_flux( + model_phot[phot_item] = syn_phot.spectrum_to_flux( model_box_full.wavelength, model_box_full.flux - ) - model_list.append(np.array([model_phot[0]])) + )[0] + model_list.append(np.array([model_phot[phot_item]])) phot_flux = object_box.flux[phot_item] phot_data = np.array( @@ -632,6 +645,7 @@ def compare_model( * model_list / data_list[:, 2] ** 2 ) + c_denom = ( weights_list * model_list**2 @@ -689,7 +703,18 @@ def compare_model( spec_idx, ] = spec_scaling - g_fit = 0.0 + for phot_item in inc_phot: + phot_flux = object_box.flux[phot_item] + + g_fit += np.nansum( + w_i[phot_item] + * ( + phot_flux[0] + - scaling * model_phot[phot_item] + ) + ** 2 + / phot_flux[1] ** 2 + ) for spec_item in self.spec_name: if spec_item in scale_spec: @@ -766,4 +791,5 @@ def compare_model( model=model, scale_spec=scale_spec, extra_scaling=extra_scaling, + inc_phot=inc_phot, ) diff --git a/species/data/database.py b/species/data/database.py index b7c4e30e..7e8723e7 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -2567,21 +2567,26 @@ def get_object( inc_spec: Union[bool, List[str]] = True, ) -> box.ObjectBox: """ - Function for extracting the photometric and/or spectroscopic data of an object from the - database. The spectroscopic data contains optionally the covariance matrix and its inverse. + Function for extracting the photometric and/or spectroscopic + data of an object from the database. The spectroscopic data + contains optionally the covariance matrix and its inverse. Parameters ---------- object_name : str Object name in the database. inc_phot : bool, list(str) - Include photometric data. If a boolean, either all (``True``) or none (``False``) of - the data are selected. If a list, a subset of filter names (as stored in the database) - can be provided. + Include photometric data. If a boolean, either all + (``True``) or none (``False``) of the data are selected. + If a list, a subset of filter names (as stored in the + database) can be provided. inc_spec : bool, list(str) - Include spectroscopic data. If a boolean, either all (``True``) or none (``False``) of - the data are selected. If a list, a subset of spectrum names (as stored in the database - with :func:`~species.data.database.Database.add_object`) can be provided. + Include spectroscopic data. If a boolean, either all + (``True``) or none (``False``) of the data are selected. + If a list, a subset of spectrum names (as stored in the + database with + :func:`~species.data.database.Database.add_object`) can + be provided. Returns ------- @@ -3064,6 +3069,7 @@ def add_comparison( model: str, scale_spec: List[str], extra_scaling: Optional[np.ndarray], + inc_phot: List[str], ) -> None: """ Function for adding results obtained with @@ -3100,6 +3106,9 @@ def add_comparison( Array with extra scalings that have been applied to the spectra of ``scale_spec``. The argument can be set to ``None`` if no extra scalings have been applied. + inc_phot : list(str) + List with filter names of which photometric data + was included with the comparison. Returns ------- @@ -3130,6 +3139,7 @@ def add_comparison( dset.attrs["n_spec_name"] = len(spec_name) dset.attrs["n_scale_spec"] = len(scale_spec) dset.attrs["parallax"] = parallax + dset.attrs["n_inc_phot"] = len(inc_phot) for i, item in enumerate(model_param): dset.attrs[f"parameter{i}"] = item @@ -3140,6 +3150,9 @@ def add_comparison( for i, item in enumerate(scale_spec): dset.attrs[f"scale_spec{i}"] = item + for i, item in enumerate(inc_phot): + dset.attrs[f"inc_phot{i}"] = item + h5_file.create_dataset( f"results/comparison/{tag}/flux_scaling", data=flux_scaling ) diff --git a/species/data/model_data.json b/species/data/model_data.json index 555280fc..5a054756 100644 --- a/species/data/model_data.json +++ b/species/data/model_data.json @@ -138,7 +138,7 @@ "wavelength range": [0.35, 250.0], "resolution": 500, "teff range": [400, 2000], - "reference": "Charney et al. (2018)", + "reference": "Charnay et al. (2018)", "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract" }, "exo-rem-highres": { @@ -148,7 +148,7 @@ "wavelength range": [0.67, 250.0], "resolution": 20000, "teff range": [400, 2000], - "reference": "Charney et al. (2018)", + "reference": "Charnay et al. (2018)", "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract" }, "koester-wd": { diff --git a/species/plot/plot_comparison.py b/species/plot/plot_comparison.py index 47ac435f..d6abb8fb 100644 --- a/species/plot/plot_comparison.py +++ b/species/plot/plot_comparison.py @@ -963,10 +963,10 @@ def plot_model_spectra( flux_offset : float, None Offset to be applied such that the spectra do not overlap. No offset is applied if the argument is set to ``None``. - xlim : tuple(float, float) - Limits of the spectral type axis. - ylim : tuple(float, float) - Limits of the goodness-of-fit axis. + xlim : tuple(float, float), None + Limits of the wavelength axis. + ylim : tuple(float, float), None + Limits of the flux axis. title : str Plot title. offset : tuple(float, float) @@ -1016,11 +1016,16 @@ def plot_model_spectra( n_param = dset.attrs["n_param"] n_scale_spec = dset.attrs["n_scale_spec"] parallax = dset.attrs["parallax"] + n_inc_phot = dset.attrs["n_inc_phot"] spec_name = [] for i in range(n_spec_name): spec_name.append(dset.attrs[f"spec_name{i}"]) + inc_phot = [] + for i in range(n_inc_phot): + inc_phot.append(dset.attrs[f"inc_phot{i}"]) + model_param = [] coord_points = [] for i in range(n_param): @@ -1174,6 +1179,19 @@ def plot_model_spectra( color="black", ) + if n_inc_phot > 0 and n_spec_name == 0: + model_box = model_reader.get_data( + model_param=param_select, + spec_res=100., + ) + + ax.plot( + model_box.wavelength, + (n_spectra - i - 1) * flux_offset + model_box.flux, + color="tomato", + lw=0.5, + ) + for spec_key, spec_item in obj_spec.items(): if spec_key not in spec_name: continue diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index 9fbea8e9..96696d84 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -1057,7 +1057,11 @@ def plot_spectrum( ) finite = np.isfinite(residuals.photometry[item][1]) - res_max = np.max(np.abs(residuals.photometry[item][1][finite])) + + max_tmp = np.max(np.abs(residuals.photometry[item][1][finite])) + + if max_tmp > res_max: + res_max = max_tmp if residuals.spectrum is not None: for key, value in residuals.spectrum.items(): diff --git a/species/util/plot_util.py b/species/util/plot_util.py index 95fd835d..c886ab96 100644 --- a/species/util/plot_util.py +++ b/species/util/plot_util.py @@ -1392,6 +1392,7 @@ def create_model_label( if model_name in model_list: read_mod = read_model.ReadModel(model_name) check_param = read_mod.get_parameters() + check_param += ['radius', 'mass', 'luminosity'] else: check_param = None