From c95ac2c57cfd2849ec4b1df0c59592ab7db46c11 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 18 Sep 2024 14:59:49 +0200 Subject: [PATCH] Several updates, functionality improvements, and fixed a few issues --- docs/tutorials/read_isochrone.ipynb | 10 +- species/core/box.py | 10 +- .../data/companion_data/companion_data.json | 17 +- species/data/database.py | 92 +-- species/fit/compare_spectra.py | 7 +- species/fit/emission_line.py | 4 +- species/fit/fit_evolution.py | 557 ++++++++++-------- species/fit/fit_model.py | 38 +- species/fit/fit_spectrum.py | 4 +- species/phot/syn_phot.py | 7 +- species/plot/plot_evolution.py | 165 ++++-- species/plot/plot_mcmc.py | 37 +- species/plot/plot_spectrum.py | 37 +- species/read/read_isochrone.py | 19 +- species/read/read_model.py | 34 +- species/util/fit_util.py | 19 +- 16 files changed, 632 insertions(+), 425 deletions(-) diff --git a/docs/tutorials/read_isochrone.ipynb b/docs/tutorials/read_isochrone.ipynb index bda8a37a..c2bef31f 100644 --- a/docs/tutorials/read_isochrone.ipynb +++ b/docs/tutorials/read_isochrone.ipynb @@ -13,7 +13,7 @@ "id": "7f6f8c20", "metadata": {}, "source": [ - "In this tutorial, we will work with data from an evolutionary model. We will extract an isochrone and cooling curve by interpolating the data at a fixed age and mass, respectively. We will also compute synthetic photometry for a given age and mass by using the associated grid with model spectra." + "In this tutorial, we will work with data from an evolutionary model. We will extract an isochrone and cooling track by interpolating the data at a fixed age and mass, respectively. We will also compute synthetic photometry for a given age and mass by using the associated grid with model spectra." ] }, { @@ -493,7 +493,7 @@ "id": "65c4ab2e", "metadata": {}, "source": [ - "## Extracting a cooling curve" + "## Extracting a cooling track" ] }, { @@ -501,7 +501,7 @@ "id": "7c6d9faa", "metadata": {}, "source": [ - "Similarly to extracting an isochrone, we can extract a cooling curve with the [get_cooling_curve](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_curve) method, which will interpolate the evolutionary data at a fixed mass and a range of ages. Instead of providing an `numpy` array with ages, we can also set the argument of `ages` to `None`. In that case it use the ages that are available in the original data." + "Similarly to extracting an isochrone, we can extract a cooling track with the [get_cooling_track](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_track) method, which will interpolate the evolutionary data at a fixed mass and a range of ages. Instead of providing an `numpy` array with ages, we can also set the argument of `ages` to `None`. In that case it use the ages that are available in the original data." ] }, { @@ -511,7 +511,7 @@ "metadata": {}, "outputs": [], "source": [ - "cooling_box = read_iso.get_cooling_curve(mass=10.,\n", + "cooling_box = read_iso.get_cooling_track(mass=10.,\n", " ages=None,\n", " filters_color=None,\n", " filter_mag=None)" @@ -522,7 +522,7 @@ "id": "a08d4f32", "metadata": {}, "source": [ - "The cooling curve that is returned by [get_cooling_curve](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_curve) is stored in a [CoolingBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.CoolingBox). Let's have a look at the content by again using the [open_box](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.Box.open_box) method." + "The cooling track that is returned by [get_cooling_track](https://species.readthedocs.io/en/latest/species.read.html#species.read.read_isochrone.ReadIsochrone.get_cooling_track) is stored in a [CoolingBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.CoolingBox). Let's have a look at the content by again using the [open_box](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.Box.open_box) method." ] }, { diff --git a/species/core/box.py b/species/core/box.py index 089f4c93..e012fb1a 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -546,13 +546,19 @@ def create_box(boxtype, **kwargs): box = ResidualsBox() box.name = kwargs["name"] box.photometry = kwargs["photometry"] - box.spectrum = kwargs["spectrum"] + if "spectrum" in kwargs: + box.spectrum = kwargs["spectrum"] + if "model_name" in kwargs: + box.model_name = kwargs["model_name"] if "chi2_red" in kwargs: box.chi2_red = kwargs["chi2_red"] elif boxtype == "samples": box = SamplesBox() - box.spectrum = kwargs["spectrum"] + if "spectrum" in kwargs: + box.spectrum = kwargs["spectrum"] + if "model_name" in kwargs: + box.model_name = kwargs["model_name"] box.parameters = kwargs["parameters"] box.samples = kwargs["samples"] box.ln_prob = kwargs["ln_prob"] diff --git a/species/data/companion_data/companion_data.json b/species/data/companion_data/companion_data.json index 47d49df0..44f3253e 100644 --- a/species/data/companion_data/companion_data.json +++ b/species/data/companion_data/companion_data.json @@ -112,7 +112,7 @@ "Currie et al. 2013", "Gaia Early Data Release 3", "Males et al. 2014", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Miret-Roig et al. 2020", "Nowak et al. 2020", "Stolker et al. 2019", @@ -152,7 +152,7 @@ ], "references": [ "Gaia Early Data Release 3", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Miret-Roig et al. 2020", "Nowak et al. 2020" ] @@ -369,7 +369,7 @@ "Gaia Early Data Release 3", "Galicher et al. 2011", "Marois et al. 2010", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Wang et al. 2018", "Zurlo et al. 2016" ] @@ -448,7 +448,7 @@ "Gaia Early Data Release 3", "Galicher et al. 2011", "Marois et al. 2010", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Wang et al. 2018", "Zurlo et al. 2016" ] @@ -527,7 +527,7 @@ "Gaia Early Data Release 3", "Galicher et al. 2011", "Marois et al. 2010", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Wang et al. 2018", "Zurlo et al. 2016" ] @@ -596,7 +596,7 @@ "Currie et al. 2014", "Gaia Early Data Release 3", "Marois et al. 2010", - "Mirek Brandt et al. 2021", + "Brandt et al. 2021", "Wang et al. 2018", "Zurlo et al. 2016" ] @@ -947,8 +947,9 @@ 0.03 ], "mass_companion": [ - 20.0, - 20.0 + 12.7, + -1.0, + 1.2 ], "accretion": false, "age": [ diff --git a/species/data/database.py b/species/data/database.py index e1a2bf52..5ab86625 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -2214,11 +2214,15 @@ def add_samples( group_path = f"results/fit/{tag}/fixed_param/{key}" hdf5_file.create_dataset(group_path, data=value) - if "spec_type" in attr_dict: - dset.attrs["type"] = attr_dict["spec_type"] + if "model_type" in attr_dict: + dset.attrs["model_type"] = attr_dict["model_type"] + elif "spec_type" in attr_dict: + dset.attrs["model_type"] = attr_dict["spec_type"] - if "spec_name" in attr_dict: - dset.attrs["spectrum"] = attr_dict["spec_name"] + if "model_name" in attr_dict: + dset.attrs["model_name"] = attr_dict["model_name"] + elif "spec_name" in attr_dict: + dset.attrs["model_name"] = attr_dict["spec_name"] dset.attrs["n_param"] = int(len(modelpar)) dset.attrs["sampler"] = str(sampler) @@ -2590,8 +2594,15 @@ def get_mcmc_spectra( hdf5_file = h5py.File(self.database, "r") dset = hdf5_file[f"results/fit/{tag}/samples"] - spectrum_type = dset.attrs["type"] - spectrum_name = dset.attrs["spectrum"] + if "model_type" in dset.attrs: + model_type = dset.attrs["model_type"] + else: + model_type = dset.attrs["type"] + + if "model_name" in dset.attrs: + model_name = dset.attrs["model_name"] + else: + model_name = dset.attrs["spectrum"] if "n_param" in dset.attrs: n_param = dset.attrs["n_param"] @@ -2635,7 +2646,7 @@ def get_mcmc_spectra( elif dset.attrs[f"parameter{i}"][:9] == "corr_amp_": ignore_param.append(dset.attrs[f"parameter{i}"]) - if spec_res is not None and spectrum_type == "calibration": + if spec_res is not None and model_type == "calibration": warnings.warn( "Smoothing of the spectral resolution is not " "implemented for calibration spectra." @@ -2678,24 +2689,24 @@ def get_mcmc_spectra( hdf5_file.close() - if spectrum_type == "model": - if spectrum_name == "planck": + if model_type in ["model", "atmosphere"]: + if model_name == "planck": from species.read.read_planck import ReadPlanck readmodel = ReadPlanck(wavel_range) - elif spectrum_name == "powerlaw": + elif model_name == "powerlaw": pass else: from species.read.read_model import ReadModel - readmodel = ReadModel(spectrum_name, wavel_range=wavel_range) + readmodel = ReadModel(model_name, wavel_range=wavel_range) - elif spectrum_type == "calibration": + elif model_type == "calibration": from species.read.read_calibration import ReadCalibration - readcalib = ReadCalibration(spectrum_name, filter_name=None) + readcalib = ReadCalibration(model_name, filter_name=None) boxes = [] @@ -2724,15 +2735,15 @@ def get_mcmc_spectra( elif "distance" not in model_param and distance is not None: model_param["distance"] = distance - if spectrum_type == "model": - if spectrum_name == "planck": + if model_type in ["model", "atmosphere"]: + if model_name == "planck": specbox = readmodel.get_spectrum( model_param, spec_res, wavel_resample=wavel_resample, ) - elif spectrum_name == "powerlaw": + elif model_name == "powerlaw": if wavel_resample is not None: warnings.warn( "The 'wavel_resample' parameter is not support by the " @@ -2779,7 +2790,7 @@ def get_mcmc_spectra( specbox = create_box( boxtype="model", - model=spectrum_name, + model=model_name, wavelength=specbox_0.wavelength, flux=flux_comb, parameters=model_param, @@ -2794,7 +2805,7 @@ def get_mcmc_spectra( ext_filter=ext_filter, ) - elif spectrum_type == "calibration": + elif model_type == "calibration": specbox = readcalib.get_spectrum(model_param) boxes.append(specbox) @@ -2864,8 +2875,15 @@ def get_mcmc_photometry( elif "nparam" in dset.attrs: n_param = dset.attrs["nparam"] - spectrum_type = dset.attrs["type"] - spectrum_name = dset.attrs["spectrum"] + if "model_type" in dset.attrs: + model_type = dset.attrs["model_type"] + else: + model_type = dset.attrs["type"] + + if "model_name" in dset.attrs: + model_name = dset.attrs["model_name"] + else: + model_name = dset.attrs["spectrum"] if "binary" in dset.attrs: binary = dset.attrs["binary"] @@ -2901,7 +2919,7 @@ def get_mcmc_photometry( for i in range(n_param): param.append(dset.attrs[f"parameter{i}"]) - if spectrum_type == "model": + if model_type in ["model", "atmosphere"]: if spectrum_name == "powerlaw": from species.phot.syn_phot import SyntheticPhotometry @@ -2913,7 +2931,7 @@ def get_mcmc_photometry( readmodel = ReadModel(spectrum_name, filter_name=filter_name) - elif spectrum_type == "calibration": + elif model_type == "calibration": from species.read.read_calibration import ReadCalibration readcalib = ReadCalibration(spectrum_name, filter_name=filter_name) @@ -2936,7 +2954,7 @@ def get_mcmc_photometry( elif "distance" not in model_param and distance is not None: model_param["distance"] = distance - if spectrum_type == "model": + if model_type in ["model", "atmosphere"]: if spectrum_name == "powerlaw": from species.util.model_util import powerlaw_spectrum @@ -3004,7 +3022,7 @@ def get_mcmc_photometry( else: mcmc_phot[i], _ = readmodel.get_flux(model_param) - elif spectrum_type == "calibration": + elif model_type == "calibration": if phot_type == "magnitude": app_mag, _ = readcalib.get_magnitude(model_param=model_param) mcmc_phot[i] = app_mag[0] @@ -3229,7 +3247,10 @@ def get_samples( for item in dset.attrs: attributes[item] = dset.attrs[item] - spectrum = dset.attrs["spectrum"] + if "model_name" in dset.attrs: + model_name = dset.attrs["model_name"] + else: + model_name = dset.attrs["spectrum"] if "n_param" in dset.attrs: n_param = dset.attrs["n_param"] @@ -3239,7 +3260,6 @@ def get_samples( if "ln_evidence" in dset.attrs: ln_evidence = dset.attrs["ln_evidence"] else: - # For backward compatibility ln_evidence = None param = [] @@ -3249,8 +3269,6 @@ def get_samples( print(f" - {param[-1]}") # Printing uniform and normal priors - # Check if attributes are present for - # backward compatibility uniform_priors = {} normal_priors = {} @@ -3337,7 +3355,7 @@ def get_samples( return create_box( "samples", - spectrum=spectrum, + model_name=model_name, parameters=param, samples=samples, ln_prob=ln_prob, @@ -3377,7 +3395,6 @@ def get_evidence(self, tag: str) -> Tuple[float, float]: if "ln_evidence" in dset.attrs: ln_evidence = dset.attrs["ln_evidence"] else: - # For backward compatibility ln_evidence = (None, None) return ln_evidence[0], ln_evidence[1] @@ -3427,12 +3444,16 @@ def get_pt_profiles( hdf5_file = h5py.File(self.database, "r") dset = hdf5_file[f"results/fit/{tag}/samples"] - spectrum = dset.attrs["spectrum"] + if "model_name" in dset.attrs: + model_name = dset.attrs["model_name"] + else: + model_name = dset.attrs["spectrum"] + pt_profile = dset.attrs["pt_profile"] - if spectrum != "petitradtrans": + if model_name != "petitradtrans": raise ValueError( - f"The model spectrum of the posterior samples is '{spectrum}' " + f"The model spectrum of the posterior samples is '{model_name}' " f"instead of 'petitradtrans'. Extracting P-T profiles is " f"therefore not possible." ) @@ -3866,8 +3887,8 @@ def add_retrieval( dset = hdf5_file.create_dataset(f"results/fit/{tag}/samples", data=samples) - dset.attrs["type"] = "model" - dset.attrs["spectrum"] = "petitradtrans" + dset.attrs["model_type"] = "retrieval" + dset.attrs["model_name"] = "petitradtrans" dset.attrs["n_param"] = len(parameters) if "parallax" in radtrans: @@ -4356,7 +4377,6 @@ def get_retrieval_spectra( temp_nodes = dset.attrs["temp_nodes"] else: - # For backward compatibility temp_nodes = None # Get distance diff --git a/species/fit/compare_spectra.py b/species/fit/compare_spectra.py index 48e75513..2debd5ca 100644 --- a/species/fit/compare_spectra.py +++ b/species/fit/compare_spectra.py @@ -268,11 +268,14 @@ def spectral_type( idx_select = np.isfinite(c_numer) & np.isfinite(c_denom) - c_k = np.sum(c_numer[idx_select]) / np.sum(c_denom[idx_select]) + c_k = np.sum(c_numer[idx_select]) / np.sum( + c_denom[idx_select] + ) c_k_spec.append(c_k) chi_sq = ( - spec_item[indices, 1][idx_select] - c_k * flux_resample[idx_select] + spec_item[indices, 1][idx_select] + - c_k * flux_resample[idx_select] ) / spec_item[indices, 2][idx_select] g_k += np.sum(w_i * chi_sq**2) diff --git a/species/fit/emission_line.py b/species/fit/emission_line.py index 62d00b9d..bc0c94ad 100644 --- a/species/fit/emission_line.py +++ b/species/fit/emission_line.py @@ -1163,8 +1163,8 @@ def lnlike_ultranest(params: np.ndarray) -> np.float64: # Dictionary with attributes that will be stored attr_dict = { - "spec_type": "model", - "spec_name": "gaussian", + "model_type": "emission line", + "model_name": "gaussian", "ln_evidence": (ln_z, ln_z_error), "parallax": self.parallax, } diff --git a/species/fit/fit_evolution.py b/species/fit/fit_evolution.py index 61d1d80e..f3c6f8ed 100644 --- a/species/fit/fit_evolution.py +++ b/species/fit/fit_evolution.py @@ -29,6 +29,7 @@ from species.core import constants from species.read.read_isochrone import ReadIsochrone +from species.util.core_util import print_section class FitEvolution: @@ -44,11 +45,12 @@ class FitEvolution: def __init__( self, evolution_model: str, - object_lbol: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]], - object_mass: Optional[ + log_lum: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]], + age_prior: Optional[Tuple[float, float, float]] = None, + mass_prior: Optional[ Union[Tuple[float, float], List[Optional[Tuple[float, float]]]] ] = None, - object_radius: Optional[ + radius_prior: Optional[ Union[Tuple[float, float], List[Optional[Tuple[float, float]]]] ] = None, bounds: Optional[ @@ -71,23 +73,28 @@ def __init__( argument, and error message is printed that includes a list with the isochrone models that are available in the current ``species`` database. - object_lbol : tuple(float, float), list(tuple(float, float)) + log_lum : tuple(float, float), list(tuple(float, float)) List with tuples that contain :math:`\\log10{L/L_\\odot}` and the related uncertainty for one or multiple objects. The list should follow the alphabetical order of companion characters (i.e. b, c, d, etc.) to make sure that the labels are correctly shown when plotting results. - object_mass : tuple(float, float), list(tuple(float, float)), None + age_prior : tuple(float, float, float), None + Tuple with an optional (asymmetric) normal prior for the + age (Myr). The tuple should contain three values, for + example, ``age_prior=(20., -5., +2.)``. The prior is not + applied if the argument is set to ``None``. + mass_prior : tuple(float, float), list(tuple(float, float)), None Optional list with tuples that contain the (dynamical) masses and the related uncertainty for one or multiple - objects. These masses we be used as Gaussian prior with - the fit. The order should be identical to ``object_lbol``. - object_radius : tuple(float, float), list(tuple(float, float)), None + objects. These masses we be used as normal prior with + the fit. The order should be identical to ``log_lum``. + radius_prior : tuple(float, float), list(tuple(float, float)), None Optional list with tuples that contain the radii (e.g. from and SED fit) and the related uncertainty for one - or multiple objects. These radii we be used as Gaussian + or multiple objects. These radii we be used as normal prior with the fit. The order should be identical to - ``object_lbol``. + ``log_lum``. bounds : dict(str, tuple(float, float)), None The boundaries that are used for the uniform or log-uniform priors. Fixing a parameter is possible by @@ -100,33 +107,44 @@ def __init__( None """ + print_section("Fit evolutionary model") + self.evolution_model = evolution_model - self.object_lbol = object_lbol - self.object_mass = object_mass - self.object_radius = object_radius + self.log_lum = log_lum + self.mass_prior = mass_prior + self.radius_prior = radius_prior self.bounds = bounds + self.age_prior = age_prior + self.normal_prior = {} self.fix_param = {} - if isinstance(self.object_lbol, tuple): - self.object_lbol = [self.object_lbol] - self.object_mass = [self.object_mass] - self.n_planets = 1 + print(f"Evolution model: {self.evolution_model}") + print(f"Luminosity log(L/Lsun): {self.log_lum}") - if isinstance(self.object_lbol, list): - self.n_planets = len(self.object_lbol) + print(f"\nAge prior: {self.age_prior}") + print(f"Mass prior (Rjup): {self.mass_prior}") + print(f"Radius prior (Rjup): {self.radius_prior}") - else: - self.n_planets = 1 + if isinstance(self.log_lum, tuple): + self.log_lum = [self.log_lum] + + self.n_planets = len(self.log_lum) - if self.object_mass is None: - self.object_mass = [] + if self.mass_prior is None: + self.mass_prior = [] for i in range(self.n_planets): - self.object_mass.append(None) + self.mass_prior.append(None) - if self.object_radius is None: - self.object_radius = [] + elif isinstance(self.mass_prior, tuple): + self.mass_prior = [self.mass_prior] + + if self.radius_prior is None: + self.radius_prior = [] for i in range(self.n_planets): - self.object_radius.append(None) + self.radius_prior.append(None) + + elif isinstance(self.radius_prior, tuple): + self.radius_prior = [self.radius_prior] config_file = os.path.join(os.getcwd(), "species_config.ini") @@ -138,9 +156,15 @@ def __init__( # Add grid with evolution data - with h5py.File(self.database_path, "r") as h5_file: - found_model = bool(f"isochrones/{evolution_model}" in h5_file) - tag_list = list(h5_file["isochrones"]) + with h5py.File(self.database_path, "r") as hdf_file: + found_group = bool("isochrones" in hdf_file) + + if found_group: + found_model = bool(f"isochrones/{evolution_model}" in hdf_file) + tag_list = list(hdf_file["isochrones"]) + else: + found_model = False + tag_list = None if not found_model: raise ValueError( @@ -159,70 +183,21 @@ def __init__( for i in range(self.n_planets): self.model_par.append(f"mass_{i}") - @typechecked - def run_multinest( - self, - tag: str, - n_live_points: int = 1000, - output: str = "multinest/", - ) -> None: - """ - Function to run the ``PyMultiNest`` wrapper of the - ``MultiNest`` sampler. While ``PyMultiNest`` can be - installed with ``pip`` from the PyPI repository, - ``MultiNest`` has to to be build manually. See the - ``PyMultiNest`` documentation for details: - http://johannesbuchner.github.io/PyMultiNest/install.html. - Note that the library path of ``MultiNest`` should be set - to the environmental variable ``LD_LIBRARY_PATH`` on a - Linux machine and ``DYLD_LIBRARY_PATH`` on a Mac. - Alternatively, the variable can be set before importing - the ``species`` package, for example: - - .. code-block:: python - - >>> import os - >>> os.environ['DYLD_LIBRARY_PATH'] = '/path/to/MultiNest/lib' - >>> import species + # Read isochrone grid - Parameters - ---------- - tag : str - Database tag where the samples will be stored. - n_live_points : int - Number of live points. - output : str - Path that is used for the output files from MultiNest. - - Returns - ------- - NoneType - None - """ - - try: - from mpi4py import MPI - - mpi_rank = MPI.COMM_WORLD.Get_rank() - - except ModuleNotFoundError: - mpi_rank = 0 - - # Create the output folder if required - - if mpi_rank == 0 and not os.path.exists(output): - os.mkdir(output) - - read_iso = ReadIsochrone( + self.read_iso = ReadIsochrone( tag=self.evolution_model, create_regular_grid=False, + verbose=False, ) # Prior boundaries if self.bounds is not None: + # Set manual prior boundaries + bounds_grid = {} - for key, value in read_iso.grid_points().items(): + for key, value in self.read_iso.grid_points().items(): if key in ["age", "mass"]: bounds_grid[key] = (value[0], value[-1]) @@ -305,9 +280,9 @@ def run_multinest( for i in range(self.n_planets): if f"inflate_mass{i}" in self.bounds: - if self.object_mass[i] is None: + if self.mass_prior[i] is None: warnings.warn( - f"The object_mass with index " + f"The mass_prior with index " f"{i} is set to None so the " f"inflate_mass{i} parameter " f"will be excluded." @@ -319,9 +294,10 @@ def run_multinest( self.model_par.append(f"inflate_mass{i}") else: - # Set all parameter boundaries to the grid boundaries + # Set all prior boundaries to the grid boundaries + self.bounds = {} - for key, value in read_iso.grid_points().items(): + for key, value in self.read_iso.grid_points().items(): if key == "age": self.bounds[key] = (value[0], value[-1]) @@ -343,20 +319,95 @@ def run_multinest( self.model_par.remove(item) if self.fix_param: - print(f"Fixing {len(self.fix_param)} parameters:") + print(f"\nFixing {len(self.fix_param)} parameters:") for key, value in self.fix_param.items(): print(f" - {key} = {value}") - print(f"Fitting {len(self.model_par)} parameters:") + print(f"\nFitting {len(self.model_par)} parameters:") for item in self.model_par: print(f" - {item}") - print("Prior boundaries:") + # Printing uniform and normal priors - for key, value in self.bounds.items(): - print(f" - {key} = {value}") + print("\nUniform priors (min, max):") + + for param_key, param_value in self.bounds.items(): + print(f" - {param_key} = {param_value}") + + if len(self.normal_prior) > 0: + print("\nNormal priors (mean, sigma):") + for param_key, param_value in self.normal_prior.items(): + if -0.1 < param_value[0] < 0.1: + print( + f" - {param_key} = {param_value[0]:.2e} +/- {param_value[1]:.2e}" + ) + else: + print( + f" - {param_key} = {param_value[0]:.2f} +/- {param_value[1]:.2f}" + ) + + @typechecked + def run_multinest( + self, + tag: str, + n_live_points: int = 200, + output: str = "multinest/", + ) -> None: + """ + Function to run the ``PyMultiNest`` wrapper of the + ``MultiNest`` sampler. While ``PyMultiNest`` can be + installed with ``pip`` from the PyPI repository, + ``MultiNest`` has to to be build manually. See the + ``PyMultiNest`` documentation for details: + http://johannesbuchner.github.io/PyMultiNest/install.html. + Note that the library path of ``MultiNest`` should be set + to the environmental variable ``LD_LIBRARY_PATH`` on a + Linux machine and ``DYLD_LIBRARY_PATH`` on a Mac. + Alternatively, the variable can be set before importing + the ``species`` package, for example: + + .. code-block:: python + + >>> import os + >>> os.environ['DYLD_LIBRARY_PATH'] = '/path/to/MultiNest/lib' + >>> import species + + Parameters + ---------- + tag : str + Database tag where the samples will be stored. + n_live_points : int + Number of live points. + output : str + Path that is used for the output files from MultiNest. + + Returns + ------- + NoneType + None + """ + + print_section("Nested sampling with MultiNest") + + print(f"Database tag: {tag}") + print(f"Number of live points: {n_live_points}") + print(f"Output folder: {output}") + print() + + try: + from mpi4py import MPI + + mpi_rank = MPI.COMM_WORLD.Get_rank() + + except ModuleNotFoundError: + mpi_rank = 0 + + # Create the output folder if required + + if mpi_rank == 0 and not os.path.exists(output): + os.mkdir(output) # Create a dictionary with the cube indices of the parameters @@ -390,9 +441,9 @@ def ln_prior(cube, n_dim, n_param) -> None: if key != "age": obj_idx = int(key.split("_")[-1]) - if key[:4] == "mass" and self.object_mass[obj_idx] is not None: - # Gaussian mass prior - sigma = self.object_mass[obj_idx][1] + if key[:4] == "mass" and self.mass_prior[obj_idx] is not None: + # Normal mass prior + sigma = self.mass_prior[obj_idx][1] if f"inflate_mass{obj_idx}" in self.bounds: sigma += ( @@ -406,7 +457,7 @@ def ln_prior(cube, n_dim, n_param) -> None: cube[cube_index[f"mass_{obj_idx}"]] = stats.norm.ppf( cube[cube_index[f"mass_{obj_idx}"]], - loc=self.object_mass[obj_idx][0], + loc=self.mass_prior[obj_idx][0], scale=sigma, ) @@ -417,7 +468,7 @@ def ln_prior(cube, n_dim, n_param) -> None: ) @typechecked - def ln_like(params, n_dim, n_param) -> np.float64: + def ln_like(params, n_dim, n_param) -> Union[float, np.float64]: """ Function for return the log-likelihood for the sampled parameter cube. @@ -441,69 +492,85 @@ def ln_like(params, n_dim, n_param) -> np.float64: chi_square = 0.0 - for i in range(self.n_planets): + for planet_idx in range(self.n_planets): param_names = [ "age", - f"mass_{i}", + f"mass_{planet_idx}", ] + if "age" in self.fix_param: + age_param = self.fix_param["age"] + else: + age_param = params[cube_index["age"]] + param_val = [] - for item in param_names: - if item in self.fix_param: - param_val.append(self.fix_param[item]) + for param_item in param_names: + if param_item in self.fix_param: + param_val.append(self.fix_param[param_item]) else: - param_val.append(params[cube_index[item]]) + param_val.append(params[cube_index[param_item]]) - if f"inflate_lbol{i}" in self.bounds: + if f"inflate_lbol{planet_idx}" in self.bounds: lbol_var = ( - self.object_lbol[i][1] + params[cube_index[f"inflate_lbol{i}"]] + self.log_lum[planet_idx][1] + + params[cube_index[f"inflate_lbol{planet_idx}"]] ) ** 2.0 else: - lbol_var = self.object_lbol[i][1] ** 2 + lbol_var = self.log_lum[planet_idx][1] ** 2 - iso_box = read_iso.get_isochrone( - age=params[cube_index["age"]], - masses=np.array([params[cube_index[f"mass_{i}"]]]), + iso_box = self.read_iso.get_isochrone( + age=age_param, + masses=np.array([params[cube_index[f"mass_{planet_idx}"]]]), filters_color=None, filter_mag=None, ) chi_square += ( - self.object_lbol[i][0] - iso_box.log_lum[0] + self.log_lum[planet_idx][0] - iso_box.log_lum[0] ) ** 2 / lbol_var # Only required when fitting the # inflation on the Lbol variance chi_square += np.log(2.0 * np.pi * lbol_var) - if self.object_mass[i] is not None: + if self.mass_prior[planet_idx] is not None: # Only required when fitting the # inflation on the mass variance - if f"inflate_mass{i}" in self.bounds: + if f"inflate_mass{planet_idx}" in self.bounds: mass_var = ( - self.object_mass[i][1] - + params[cube_index[f"inflate_mass{i}"]] + self.mass_prior[planet_idx][1] + + params[cube_index[f"inflate_mass{planet_idx}"]] ) ** 2.0 else: - mass_var = self.object_mass[i][1] ** 2 + mass_var = self.mass_prior[planet_idx][1] ** 2 chi_square += np.log(2.0 * np.pi * mass_var) # Radius prior - if self.object_radius[i] is not None: + if self.radius_prior[planet_idx] is not None: chi_square += ( - self.object_radius[i][0] - iso_box.radius[0] - ) ** 2 / self.object_radius[i][1] ** 2 + self.radius_prior[planet_idx][0] - iso_box.radius[0] + ) ** 2 / self.radius_prior[planet_idx][1] ** 2 + + # Age prior + if self.age_prior is not None: + if age_param < self.age_prior[0]: + # Use lower errorbar on the age + chi_square += ( + self.age_prior[0] - age_param + ) ** 2 / self.age_prior[1] ** 2 + else: + # Use upper errorbar on the age + chi_square += ( + self.age_prior[0] - age_param + ) ** 2 / self.age_prior[2] ** 2 # ln_like += -0.5 * weight * (obj_item[0] - phot_flux) ** 2 / phot_var # ln_like += -0.5 * weight * np.log(2.0 * np.pi * phot_var) - if np.isnan(chi_square): - log_like = -np.inf - - elif np.isinf(chi_square): + if not np.isfinite(chi_square): log_like = -np.inf else: @@ -533,9 +600,9 @@ def ln_like(params, n_dim, n_param) -> np.float64: # Nested sampling global log-evidence ln_z = sampling_stats["nested sampling global log-evidence"] ln_z_error = sampling_stats["nested sampling global log-evidence error"] - print(f"Nested sampling global log-evidence: {ln_z:.2f} +/- {ln_z_error:.2f}") + print(f"\nNested sampling global log-evidence: {ln_z:.2f} +/- {ln_z_error:.2f}") - # Nested sampling global log-evidence + # Nested importance sampling global log-evidence ln_z = sampling_stats["nested importance sampling global log-evidence"] ln_z_error = sampling_stats[ "nested importance sampling global log-evidence error" @@ -544,127 +611,106 @@ def ln_like(params, n_dim, n_param) -> np.float64: f"Nested importance sampling global log-evidence: {ln_z:.2f} +/- {ln_z_error:.2f}" ) - # Get the best-fit (highest likelihood) point - print("Sample with the highest likelihood:") - best_params = analyzer.get_best_fit() + # Get the maximum likelihood sample + best_params = analyzer.get_best_fit() max_lnlike = best_params["log_likelihood"] + + print("\nSample with the maximum likelihood:") print(f" - Log-likelihood = {max_lnlike:.2f}") - for i, item in enumerate(best_params["parameters"]): - print(f" - {self.model_par[i]} = {item:.2f}") + for param_idx, param_item in enumerate(best_params["parameters"]): + if -0.1 < param_item < 0.1: + print(f" - {self.model_par[param_idx]} = {param_item:.2e}") + else: + print(f" - {self.model_par[param_idx]} = {param_item:.2f}") # Get the posterior samples - samples = analyzer.get_equal_weighted_posterior() - analyzer = pymultinest.analyse.Analyzer( - len(self.model_par), outputfiles_basename=output - ) + post_samples = analyzer.get_equal_weighted_posterior() - sampling_stats = analyzer.get_stats() + # Samples and ln(L) - ln_z = sampling_stats["nested sampling global log-evidence"] - ln_z_error = sampling_stats["nested sampling global log-evidence error"] - print(f"Nested sampling global log-evidence: {ln_z:.2f} +/- {ln_z_error:.2f}") - - ln_z = sampling_stats["nested importance sampling global log-evidence"] - ln_z_error = sampling_stats[ - "nested importance sampling global log-evidence error" - ] - print( - f"Nested importance sampling global log-evidence: {ln_z:.2f} +/- {ln_z_error:.2f}" - ) - - print("Sample with the highest likelihood:") - best_params = analyzer.get_best_fit() - - max_lnlike = best_params["log_likelihood"] - print(f" - Log-likelihood = {max_lnlike:.2f}") - - for i, item in enumerate(best_params["parameters"]): - print(f" - {self.model_par[i]} = {item:.2f}") - - samples = analyzer.get_equal_weighted_posterior() - - ln_prob = samples[:, -1] + ln_prob = post_samples[:, -1] + samples = post_samples[:, :-1] # Adding the fixed parameters to the samples if self.fix_param: - samples_tmp = samples[:, :-1] + samples_tmp = np.copy(samples) self.model_par = ["age"] - for i in range(self.n_planets): - self.model_par.append(f"mass_{i}") + for planet_idx in range(self.n_planets): + self.model_par.append(f"mass_{planet_idx}") - for i in range(self.n_planets): - if f"inflate_lbol{i}" in self.bounds: - self.model_par.append(f"inflate_lbol{i}") + for planet_idx in range(self.n_planets): + if f"inflate_lbol{planet_idx}" in self.bounds: + self.model_par.append(f"inflate_lbol{planet_idx}") - for i in range(self.n_planets): - if f"inflate_mass{i}" in self.bounds: - self.model_par.append(f"inflate_mass{i}") + for planet_idx in range(self.n_planets): + if f"inflate_mass{planet_idx}" in self.bounds: + self.model_par.append(f"inflate_mass{planet_idx}") - samples = np.zeros((samples_tmp.shape[0], len(self.model_par))) + samples = np.zeros((samples.shape[0], len(self.model_par))) - for i, key in enumerate(self.model_par): - if key in self.fix_param: - samples[:, i] = np.full(samples_tmp.shape[0], self.fix_param[key]) + for param_idx, param_item in enumerate(self.model_par): + if param_item in self.fix_param: + samples[:, param_idx] = np.full( + samples_tmp.shape[0], self.fix_param[param_item] + ) else: - samples[:, i] = samples_tmp[:, cube_index[key]] - - else: - samples = samples[:, :-1] + samples[:, param_idx] = samples_tmp[:, cube_index[param_item]] - # Recreate cube_index dictionary because of included fix_param + # Recreate cube_index dictionary because the fix_param + # parameters have been included in the samples array cube_index = {} - for i, item in enumerate(self.model_par): - cube_index[item] = i + for param_idx, param_item in enumerate(self.model_par): + cube_index[param_item] = param_idx - # Add atmospheric parameters (R, Teff, and log(g)) + # Add atmospheric parameters: R, Teff, and log(g) - print("Extracting the posteriors of Teff, R, and log(g)...") + print("\nExtracting the posteriors of Teff, R, and log(g):") radius = np.zeros((samples.shape[0], self.n_planets)) log_g = np.zeros((samples.shape[0], self.n_planets)) t_eff = np.zeros((samples.shape[0], self.n_planets)) - for j in tqdm(range(self.n_planets)): - for i in tqdm(range(samples.shape[0]), leave=False): - age = samples[i, cube_index["age"]] - mass = samples[i, cube_index[f"mass_{j}"]] + for planet_idx in tqdm(range(self.n_planets)): + for sample_idx in tqdm(range(samples.shape[0]), leave=False): + age = samples[sample_idx, cube_index["age"]] + mass = samples[sample_idx, cube_index[f"mass_{planet_idx}"]] - iso_box = read_iso.get_isochrone( + iso_box = self.read_iso.get_isochrone( age=age, masses=np.array([mass]), filters_color=None, filter_mag=None, ) - radius[i, j] = iso_box.radius[0] - log_g[i, j] = iso_box.logg[0] + radius[sample_idx, planet_idx] = iso_box.radius[0] + log_g[sample_idx, planet_idx] = iso_box.logg[0] l_bol = 10.0 ** iso_box.log_lum[0] * constants.L_SUN - t_eff[i, j] = ( + t_eff[sample_idx, planet_idx] = ( l_bol / ( 4.0 * np.pi - * (radius[i, j] * constants.R_JUP) ** 2 + * (radius[sample_idx, planet_idx] * constants.R_JUP) ** 2 * constants.SIGMA_SB ) ) ** 0.25 - for i in range(self.n_planets): - self.model_par.append(f"teff_evol_{i}") + for planet_idx in range(self.n_planets): + self.model_par.append(f"teff_{planet_idx}") - for i in range(self.n_planets): - self.model_par.append(f"radius_evol_{i}") + for planet_idx in range(self.n_planets): + self.model_par.append(f"radius_{planet_idx}") - for i in range(self.n_planets): - self.model_par.append(f"logg_evol_{i}") + for planet_idx in range(self.n_planets): + self.model_par.append(f"logg_{planet_idx}") samples = np.hstack((samples, t_eff, radius, log_g)) @@ -672,8 +718,8 @@ def ln_like(params, n_dim, n_param) -> np.float64: # derived parameters that were included cube_index = {} - for i, item in enumerate(self.model_par): - cube_index[item] = i + for param_idx, param_item in enumerate(self.model_par): + cube_index[param_item] = param_idx # Remove outliers @@ -697,73 +743,93 @@ def ln_like(params, n_dim, n_param) -> np.float64: # Apply uncertainty inflation - for i in range(self.n_planets): - if f"inflate_lbol{i}" in self.bounds: - # sigma_add = np.median(samples[:, cube_index[f"inflate_lbol{i}"]]) + for planet_idx in range(self.n_planets): + if f"inflate_lbol{planet_idx}" in self.bounds: + # sigma_add = np.median(samples[:, cube_index[f"inflate_lbol{planet_idx}"]]) index_prob = np.argmax(ln_prob) - sigma_add = samples[index_prob, cube_index[f"inflate_lbol{i}"]] + sigma_add = samples[index_prob, cube_index[f"inflate_lbol{planet_idx}"]] - self.object_lbol[i] = ( - self.object_lbol[i][0], - self.object_lbol[i][1] + sigma_add, + self.log_lum[planet_idx] = ( + self.log_lum[planet_idx][0], + self.log_lum[planet_idx][1] + sigma_add, ) - for i in range(self.n_planets): - if f"inflate_mass{i}" in self.bounds: - # sigma_add = np.median(samples[:, cube_index[f"inflate_lbol{i}"]]) + for planet_idx in range(self.n_planets): + if f"inflate_mass{planet_idx}" in self.bounds: + # sigma_add = np.median(samples[:, cube_index[f"inflate_lbol{planet_idx}"]]) index_prob = np.argmax(ln_prob) - sigma_add = samples[index_prob, cube_index[f"inflate_mass{i}"]] + sigma_add = samples[index_prob, cube_index[f"inflate_mass{planet_idx}"]] - self.object_mass[i] = ( - self.object_mass[i][0], - self.object_mass[i][1] + sigma_add, + self.mass_prior[planet_idx] = ( + self.mass_prior[planet_idx][0], + self.mass_prior[planet_idx][1] + sigma_add, ) - # Set object_radius to posterior value if argument was None + # Set radius_prior to posterior value if argument was None - for i, item in enumerate(self.object_radius): - if item is None: - radius_samples = samples[:, cube_index[f"radius_evol_{i}"]] - self.object_radius[i] = ( - np.mean(radius_samples), - np.std(radius_samples), - ) + # for planet_idx, planet_item in enumerate(self.radius_prior): + # if planet_item is None: + # radius_samples = samples[:, cube_index[f"radius_{planet_idx}"]] + # self.radius_prior[planet_idx] = ( + # np.mean(radius_samples), + # np.std(radius_samples), + # ) - # Adjust object_mass to posterior value + # Adjust mass_prior to posterior value - # for i, item in enumerate(self.object_mass): + # for i, item in enumerate(self.mass_prior): # mass_samples = samples[:, cube_index[f"mass_{i}"]] - # self.object_mass[i] = (np.mean(mass_samples), np.std(mass_samples)) + # self.mass_prior[i] = (np.mean(mass_samples), np.std(mass_samples)) + + # Adjust radius_prior to posterior value + + # for i, item in enumerate(self.radius_prior): + # radius_samples = samples[:, cube_index[f"radius_{i}"]] + # self.radius_prior[i] = (np.mean(radius_samples), np.std(radius_samples)) - # Adjust object_radius to posterior value + # Set age_prior to NaN if no prior was provided - # for i, item in enumerate(self.object_radius): - # radius_samples = samples[:, cube_index[f"radius_evol_{i}"]] - # self.object_radius[i] = (np.mean(radius_samples), np.std(radius_samples)) + if self.age_prior is None: + self.age_prior = [np.nan] - # Set object_mass and object_radius to NaN if no prior was provided + elif "age" in self.fix_param: + self.age_prior = (self.fix_param["age"], 0.0) + + # Set mass_prior to NaN if no prior was provided + + for planet_idx, planet_item in enumerate(self.mass_prior): + if f"mass_{planet_idx}" in self.fix_param: + self.mass_prior[planet_idx] = ( + self.fix_param[f"mass_{planet_idx}"], + 0.0, + ) - for i, item in enumerate(self.object_mass): - if f"mass_{i}" in self.fix_param: - self.object_mass[i] = (self.fix_param[f"mass_{i}"], 0.0) - elif item is None: - self.object_mass[i] = np.nan + elif planet_item is None: + self.mass_prior[planet_idx] = np.nan - for i, item in enumerate(self.object_radius): - if item is None: - self.object_radius[i] = np.nan + # Set radius_prior to NaN if no prior was provided + + for planet_idx, planet_item in enumerate(self.radius_prior): + if f"radius_{planet_idx}" in self.fix_param: + self.radius_prior[planet_idx] = ( + self.fix_param[f"radius_{planet_idx}"], + 0.0, + ) + + elif planet_item is None: + self.radius_prior[planet_idx] = np.nan # Dictionary with attributes that will be stored attr_dict = { - "spec_type": "model", - "spec_name": "evolution", - "evolution_model": self.evolution_model, + "model_type": "evolution", + "model_name": self.evolution_model, "ln_evidence": (ln_z, ln_z_error), "n_planets": self.n_planets, - "object_lbol": self.object_lbol, - "object_mass": self.object_mass, - "object_radius": self.object_radius, + "log_lum": self.log_lum, + "age_prior": self.age_prior, + "mass_prior": self.mass_prior, + "radius_prior": self.radius_prior, } # Add samples to the database @@ -776,10 +842,13 @@ def ln_like(params, n_dim, n_param) -> np.float64: species_db = Database() species_db.add_samples( + tag=tag, sampler="multinest", samples=samples, ln_prob=ln_prob, - tag=tag, modelpar=self.model_par, + bounds=self.bounds, + normal_prior=self.normal_prior, + fixed_param=self.fix_param, attr_dict=attr_dict, ) diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py index 212d383b..24b63a49 100644 --- a/species/fit/fit_model.py +++ b/species/fit/fit_model.py @@ -34,7 +34,6 @@ from scipy.stats import norm from typeguard import typechecked -from species.core import constants from species.phot.syn_phot import SyntheticPhotometry from species.read.read_model import ReadModel from species.read.read_object import ReadObject @@ -43,10 +42,8 @@ from species.util.convert_util import logg_to_mass from species.util.core_util import print_section from species.util.dust_util import ( - convert_to_av, interp_lognorm, interp_powerlaw, - ism_extinction, ) from species.util.model_util import ( binary_to_single, @@ -83,7 +80,10 @@ def __init__( str, Union[ Optional[Tuple[Optional[float], Optional[float]]], - Tuple[Optional[Tuple[Optional[float], Optional[float]]], Optional[Tuple[Optional[float], Optional[float]]]], + Tuple[ + Optional[Tuple[Optional[float], Optional[float]]], + Optional[Tuple[Optional[float], Optional[float]]], + ], Tuple[ Optional[Tuple[float, float]], Optional[Tuple[float, float]], @@ -944,9 +944,12 @@ def __init__( for spec_key, spec_value in self.spectrum.items(): print(f"\rInterpolating {spec_key}...", end="", flush=True) + # The ReadModel class will make sure that wavelength range of + # the selected subgrid fully includes the requested wavel_range + wavel_range = ( - 0.9 * spec_value[0][0, 0], - 1.1 * spec_value[0][-1, 0], + spec_value[0][0, 0], + spec_value[0][-1, 0], ) readmodel = ReadModel(self.model, wavel_range=wavel_range) @@ -1019,8 +1022,16 @@ def __init__( print(" [DONE]") for spec_key, spec_value in self.spectrum.items(): + # Use the wavelength range of the selected atmospheric model + # because the sampled blackbody spectrum will get interpolated + # to the wavelengths of the model spectrum instead of the data + spec_idx = list(self.spectrum.keys()).index(spec_key) print(f"\rInterpolating {spec_key}...", end="", flush=True) - wavel_range = (0.9 * spec_value[0][0, 0], 1.1 * spec_value[0][-1, 0]) + # wavel_range = (spec_value[0][0, 0], spec_value[0][-1, 0]) + wavel_range = ( + self.modelspec[spec_idx].wl_points[0], + self.modelspec[spec_idx].wl_points[-1], + ) readmodel = ReadModel("blackbody", wavel_range=wavel_range) readmodel.interpolate_grid(teff_range=None) self.diskspec.append(readmodel) @@ -1169,7 +1180,7 @@ def __init__( for item in self.modelpar: print(f" - {item}") - # Add parallax to dictionary with Gaussian priors + # Add parallax to dictionary with normal priors if ( "parallax" in self.modelpar @@ -1523,7 +1534,8 @@ def _lnlike_func( if phot_flux_1 / phot_flux_0 < ratio_prior[0]: return -np.inf - elif phot_flux_1 / phot_flux_0 > ratio_prior[1]: + + if phot_flux_1 / phot_flux_0 > ratio_prior[1]: return -np.inf # Normal prior for the flux ratio @@ -2057,8 +2069,8 @@ def _create_attr_dict(self): """ attr_dict = { - "spec_type": "model", - "spec_name": self.model, + "model_type": "atmosphere", + "model_name": self.model, "ln_evidence": (self.ln_z, self.ln_z_error), "parallax": self.obj_parallax[0], "binary": self.binary, @@ -2297,8 +2309,12 @@ def _lnlike_multinest( print(f" - {self.modelpar[param_idx]} = {param_item:.2f}") # Get the posterior samples + post_samples = analyzer.get_equal_weighted_posterior() + # Create list with the spectrum labels that are used + # for fitting an additional scaling parameter + spec_labels = [] for spec_item in self.spectrum: if f"scaling_{spec_item}" in self.bounds: diff --git a/species/fit/fit_spectrum.py b/species/fit/fit_spectrum.py index 479da123..b4fae113 100644 --- a/species/fit/fit_spectrum.py +++ b/species/fit/fit_spectrum.py @@ -214,8 +214,8 @@ def run_mcmc( # Dictionary with attributes that will be stored attr_dict = { - "spec_type": "calibration", - "spec_name": self.spectrum, + "model_type": "calibration", + "model_name": self.spectrum, "mean_accept": np.mean(ens_sampler.acceptance_fraction), } diff --git a/species/phot/syn_phot.py b/species/phot/syn_phot.py index 2764a247..3ce364ac 100644 --- a/species/phot/syn_phot.py +++ b/species/phot/syn_phot.py @@ -115,10 +115,7 @@ def calc_zero_point(self) -> Union[float, np.float64]: self.wavel_range = read_filt.wavelength_range() with h5py.File(self.database, "r") as hdf5_file: - if "spectra/calibration/vega" in hdf5_file: - vega_found = True - else: - vega_found = False + vega_found = "spectra/calibration/vega" in hdf5_file if not vega_found: with h5py.File(self.database, "a") as hdf5_file: @@ -285,7 +282,7 @@ def spectrum_to_flux( f"The filter profile of {self.filter_name} " f"({self.wavel_range[0]:.4f}-{self.wavel_range[1]:.4f}) " f"lies outside the wavelength range of the spectrum. " - f"The synthetic flux is set to NaN." + f"The returned synthetic flux is therefore set to NaN." ) syn_flux = np.nan diff --git a/species/plot/plot_evolution.py b/species/plot/plot_evolution.py index 25056742..74654dd9 100644 --- a/species/plot/plot_evolution.py +++ b/species/plot/plot_evolution.py @@ -13,6 +13,7 @@ from matplotlib.ticker import AutoMinorLocator from species.read.read_isochrone import ReadIsochrone +from species.util.core_util import print_section @typechecked @@ -28,9 +29,9 @@ def plot_cooling( offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, -) -> Tuple[mpl.figure.Figure, List[List[np.ndarray]], np.ndarray]: +) -> Tuple[mpl.figure.Figure, List[List[List[np.ndarray]]], np.ndarray]: """ - Function for plotting samples of cooling curves that are + Function for plotting samples of cooling tracks that are randomly drawn from the posterior distributions of the age and mass parameters that have been estimated with :class:`~species.fit.fit_evolution.FitEvolution`. @@ -40,7 +41,7 @@ def plot_cooling( tag : str Database tag where the samples are stored n_samples : int - Number of randomly drawn cooling curves that will be plotted. + Number of randomly drawn cooling tracks that will be plotted. cooling_param : str Type of cooling parameter that will be plotted ('luminosity' or 'radius'). @@ -72,7 +73,7 @@ def plot_cooling( The ``Figure`` object that can be used for further customization of the plot. np.ndarray - Array with the cooling curves. The array contains + Array with the cooling tracks. The array contains :math:`L/L_\\odot` or radius as function of time for each companion and sample, so the shape is (n_companions, n_samples, n_ages). @@ -81,6 +82,12 @@ def plot_cooling( sampled from the posterior distribution. """ + print_section("Plot cooling tracks") + + print(f"Database tag: {tag}") + print(f"Number of samples: {n_samples}") + print(f"Model parameters: {cooling_param}") + plt.close() if cooling_param not in ["luminosity", "radius"]: @@ -98,25 +105,27 @@ def plot_cooling( samples = samples_box.samples attr = samples_box.attributes n_planets = attr["n_planets"] - evolution_model = attr["evolution_model"] - object_lbol = attr["object_lbol"] - object_radius = attr["object_radius"] - object_age = (np.mean(samples[:, 0]), np.std(samples[:, 0])) + model_name = attr["model_name"] + log_lum = attr["log_lum"] + age_prior = attr["age_prior"] + radius_prior = attr["radius_prior"] + + if np.isnan(age_prior[0]): + param_idx = samples_box.parameters.index("age") + age_prior = np.percentile(samples[:, param_idx], [50.0, 16.0, 84.0]) + else: + # Asymmetric normal prior set in FitEvolution + age_prior = [ + age_prior[0], + age_prior[0] + age_prior[1], + age_prior[0] + age_prior[2], + ] - read_iso = ReadIsochrone(evolution_model) + read_iso = ReadIsochrone(model_name) plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" - if output is None: - print("Plotting cooling curves...", end="", flush=True) - else: - print( - f"Plotting cooling curves: {output}...", - end="", - flush=True, - ) - fig = plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(n_planets, 1) gridsp.update(wspace=0, hspace=0.1, left=0, right=1, bottom=0, top=1) @@ -131,11 +140,21 @@ def plot_cooling( if yscale is None: yscale = "linear" - cool_curves = [] + cool_tracks = [] for i in range(n_planets): - cool_curves.append([]) + cool_tracks.append([]) for i in range(n_planets): + if not isinstance(radius_prior[i], np.ndarray) and np.isnan(radius_prior[i]): + param_idx = samples_box.parameters.index(f"radius_{i}") + radius_tmp = np.percentile(samples[:, param_idx], [50.0, 16.0, 84.0]) + else: + radius_tmp = [ + radius_prior[i][0], + radius_prior[i][0] - radius_prior[i][1], + radius_prior[i][0] + radius_prior[i][1], + ] + ax[i].set_xscale(xscale) ax[i].set_yscale(yscale) @@ -201,17 +220,17 @@ def plot_cooling( for planet_idx in range(n_planets): mass = samples[sample_idx, 1 + planet_idx] - cool_box = read_iso.get_cooling_curve(mass=mass, ages=None) + cool_box = read_iso.get_cooling_track(mass=mass, ages=None) if cooling_param == "luminosity": - cool_curves[planet_idx].append(cool_box.log_lum) + cool_tracks[planet_idx].append([cool_box.age, cool_box.log_lum]) elif cooling_param == "radius": - cool_curves[planet_idx].append(cool_box.radius) + cool_tracks[planet_idx].append([cool_box.age, cool_box.radius]) ax[planet_idx].plot( - cool_box.age, - cool_curves[planet_idx][-1], + cool_tracks[planet_idx][-1][0], + cool_tracks[planet_idx][-1][1], lw=0.5, color="gray", alpha=0.5, @@ -220,21 +239,25 @@ def plot_cooling( for i in range(n_planets): if cooling_param == "luminosity": ax[i].errorbar( - object_age[0], - object_lbol[i][0], - xerr=object_age[1], - yerr=object_lbol[i][1], + [age_prior[0]], + [log_lum[i][0]], + xerr=[ + [age_prior[0] - np.abs(age_prior[1])], + [age_prior[2] - age_prior[0]], + ], + yerr=[log_lum[i][1]], color="tab:orange", ) - elif cooling_param == "radius" and isinstance(object_radius[i], np.ndarray): - # Only plot the data if these were provided as optional - # argument of object_radius when using FitEvolution + elif cooling_param == "radius": ax[i].errorbar( - object_age[0], - object_radius[i][0], - xerr=object_age[1], - yerr=object_radius[i][1], + [age_prior[0]], + [radius_tmp[0]], + xerr=[ + [age_prior[0] - np.abs(age_prior[1])], + [age_prior[2] - age_prior[0]], + ], + yerr=[[radius_tmp[0] - radius_tmp[1]], [radius_tmp[2] - radius_tmp[0]]], color="tab:orange", ) @@ -244,11 +267,10 @@ def plot_cooling( if output is None: plt.show() else: + print(f"\nOutput: {output}") plt.savefig(output, bbox_inches="tight") - print(" [DONE]") - - return fig, cool_curves, ran_indices + return fig, cool_tracks, ran_indices @typechecked @@ -263,7 +285,7 @@ def plot_isochrones( offset: Optional[Tuple[float, float]] = None, figsize: Optional[Tuple[float, float]] = (4.0, 2.5), output: Optional[str] = None, -) -> Tuple[mpl.figure.Figure, List[List[np.ndarray]], np.ndarray]: +) -> Tuple[mpl.figure.Figure, List[List[List[np.ndarray]]], np.ndarray]: """ Function for plotting samples of isochrones that are randomly drawn from the posterior distributions of the @@ -277,7 +299,7 @@ def plot_isochrones( tag : str Database tag where the samples are stored n_samples : int - Number of randomly drawn cooling curves that will be plotted. + Number of randomly drawn cooling tracks that will be plotted. xlim : tuple(float, float), None Limits of the wavelength axis. Automatic limits are used if the argument is set to ``None``. @@ -314,6 +336,11 @@ def plot_isochrones( sampled from the posterior distribution. """ + print_section("Plot isochrones") + + print(f"Database tag: {tag}") + print(f"Number of samples: {n_samples}") + plt.close() from species.data.database import Database @@ -324,24 +351,27 @@ def plot_isochrones( samples = samples_box.samples attr = samples_box.attributes n_planets = attr["n_planets"] - evolution_model = attr["evolution_model"] - object_lbol = attr["object_lbol"] - object_mass = attr["object_mass"] + model_name = attr["model_name"] + log_lum = attr["log_lum"] + age_prior = attr["age_prior"] + mass_prior = attr["mass_prior"] + + if np.isnan(age_prior[0]): + param_idx = samples_box.parameters.index("age") + age_prior = np.percentile(samples[:, param_idx], [50.0, 16.0, 84.0]) + else: + # Asymmetric normal prior set in FitEvolution + age_prior = [ + age_prior[0], + age_prior[0] + age_prior[1], + age_prior[0] + age_prior[2], + ] - read_iso = ReadIsochrone(evolution_model) + read_iso = ReadIsochrone(model_name) plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" - if output is None: - print("Plotting isochrones...", end="", flush=True) - else: - print( - f"Plotting isochrones: {output}...", - end="", - flush=True, - ) - fig = plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(n_planets, 1) gridsp.update(wspace=0, hspace=0.1, left=0, right=1, bottom=0, top=1) @@ -424,22 +454,32 @@ def plot_isochrones( iso_box = read_iso.get_isochrone(age=age, masses=None) - isochrones[planet_idx].append(iso_box.log_lum) + isochrones[planet_idx].append([iso_box.mass, iso_box.log_lum]) ax[planet_idx].plot( - iso_box.mass, - isochrones[planet_idx][-1], + isochrones[planet_idx][-1][0], + isochrones[planet_idx][-1][1], lw=0.5, color="gray", alpha=0.5, ) for i in range(n_planets): + if not isinstance(mass_prior[i], np.ndarray) and np.isnan(mass_prior[i]): + param_idx = samples_box.parameters.index(f"mass_{i}") + mass_tmp = np.percentile(samples[:, param_idx], [50.0, 16.0, 84.0]) + else: + mass_tmp = [ + mass_prior[i][0], + mass_prior[i][0] - mass_prior[i][1], + mass_prior[i][0] + mass_prior[i][1], + ] + ax[i].errorbar( - object_mass[i][0], - object_lbol[i][0], - xerr=object_mass[i][1], - yerr=object_lbol[i][1], + [mass_tmp[0]], + [log_lum[i][0]], + xerr=[[mass_tmp[0] - mass_tmp[1]], [mass_tmp[2] - mass_tmp[0]]], + yerr=[log_lum[i][1]], color="tab:orange", ) @@ -447,10 +487,9 @@ def plot_isochrones( ax[0].set_title(title, fontsize=18.0) if output is None: + print(f"\nOutput: {output}") plt.show() else: plt.savefig(output, bbox_inches="tight") - print(" [DONE]") - return fig, isochrones, ran_indices diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py index 87354334..3a85ef3f 100644 --- a/species/plot/plot_mcmc.py +++ b/species/plot/plot_mcmc.py @@ -313,7 +313,20 @@ def plot_posterior( print(f"Database tag: {tag}") print(f"Object type: {object_type}") - print(f"Manual parameters: {param_inc}") + print(f"Manual parameters: {param_inc}\n") + + if "model_type" in box.attributes: + print((f"Model type: {box.attributes['model_type']}")) + elif "spec_type" in box.attributes: + print((f"Model type: {box.attributes['spec_type']}")) + + if "model_name" in box.attributes: + print((f"Model name: {box.attributes['model_name']}")) + elif "spec_name" in box.attributes: + print((f"Model type: {box.attributes['spec_name']}")) + + if "sampler" in box.attributes: + print((f"Sampler: {box.attributes['sampler']}")) plt.rcParams["font.family"] = "serif" plt.rcParams["mathtext.fontset"] = "dejavuserif" @@ -971,6 +984,20 @@ def plot_posterior( ndim = len(param_inc) samples = param_new + # Only for fitting evolutionary models + # Remove index from parameter names when fitting 1 planet + + if "model_type" in box.attributes and box.attributes["model_type"] == "evolution": + if box.attributes["n_planets"] == 1: + param_copy = box.parameters.copy() + box.parameters = [] + + for param_item in param_copy: + if param_item[-2:] == "_0": + box.parameters.append(param_item[:-2]) + else: + box.parameters.append(param_item) + # Update axes labels box_param = box.parameters.copy() @@ -1000,7 +1027,7 @@ def plot_posterior( if max_prob: max_idx = np.argmax(box.ln_prob) - max_sample = samples[max_idx, ] + max_sample = samples[max_idx,] if isinstance(title_fmt, list) and len(title_fmt) != ndim: raise ValueError( @@ -1042,6 +1069,8 @@ def plot_posterior( hist_titles.append(hist_title) + # Create corner plot + fig = corner.corner( samples, quantiles=[0.16, 0.5, 0.84], @@ -1148,12 +1177,10 @@ def plot_posterior( if title: fig.suptitle(title, y=1.02, fontsize=16) - if output is not None: - print(f"\nOutput: {output}") - if output is None: plt.show() else: + print(f"\nOutput: {output}") plt.savefig(output, bbox_inches="tight") return fig diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index 4e4c95e5..0016da3a 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -425,16 +425,16 @@ def plot_spectrum( font_size = {} if "xlabel" not in font_size: - font_size["xlabel"] = 11. + font_size["xlabel"] = 11.0 if "ylabel" not in font_size: - font_size["ylabel"] = 11. + font_size["ylabel"] = 11.0 if "title" not in font_size: - font_size["title"] = 13. + font_size["title"] = 13.0 if "legend" not in font_size: - font_size["legend"] = 9. + font_size["legend"] = 9.0 print(f"Font sizes: {font_size}") @@ -459,10 +459,14 @@ def plot_spectrum( if residuals is not None: if quantity == "flux density": - ax3.set_ylabel(r"$\Delta$$F_\lambda$ ($\sigma$)", fontsize=font_size["ylabel"]) + ax3.set_ylabel( + r"$\Delta$$F_\lambda$ ($\sigma$)", fontsize=font_size["ylabel"] + ) elif quantity == "flux": - ax3.set_ylabel(r"$\Delta$$F_\lambda$ ($\sigma$)", fontsize=font_size["ylabel"]) + ax3.set_ylabel( + r"$\Delta$$F_\lambda$ ($\sigma$)", fontsize=font_size["ylabel"] + ) if quantity == "magnitude": scaling = 1.0 @@ -521,7 +525,9 @@ def plot_spectrum( ) elif quantity == "flux": - ax1.set_ylabel(r"$\lambda$$F_\lambda$ (W m$^{-2}$)", fontsize=font_size["ylabel"]) + ax1.set_ylabel( + r"$\lambda$$F_\lambda$ (W m$^{-2}$)", fontsize=font_size["ylabel"] + ) scaling = 1.0 @@ -641,7 +647,10 @@ def plot_spectrum( flux_scaling = wavelength ax1.plot( - wavelength, flux_scaling * flux_masked / scaling, lw=0.5, label=label + wavelength, + flux_scaling * flux_masked / scaling, + lw=0.5, + label=label, ) elif isinstance(box_item, list): @@ -960,7 +969,7 @@ def plot_spectrum( # wavelengths but hwhm_up and hwhm_down will # be different when converting a FWHM from # wavelength to frequency - fwhm = (hwhm_up + hwhm_down) + fwhm = hwhm_up + hwhm_down if not plot_kwargs[j] or filter_item not in plot_kwargs[j]: if not plot_kwargs[j]: @@ -1106,18 +1115,14 @@ def plot_spectrum( flux_conv = data_out[:, 1] # Convert FWHM of filter to requested units - data_in = np.column_stack( - [[wavel_micron + fwhm_micron / 2.0], [1.0]] - ) + data_in = np.column_stack([[wavel_micron + fwhm_micron / 2.0], [1.0]]) data_out = convert_units(data_in, units, convert_from=False) # Absolute value because could be negative when frequency hwhm_up = np.abs(data_out[0, 0] - wavelength[0]) # Convert FWHM of filter to requested units - data_in = np.column_stack( - [[wavel_micron - fwhm_micron / 2.0], [1.0]] - ) + data_in = np.column_stack([[wavel_micron - fwhm_micron / 2.0], [1.0]]) data_out = convert_units(data_in, units, convert_from=False) # Absolute value because could be negative when frequency @@ -1128,7 +1133,7 @@ def plot_spectrum( # wavelengths but hwhm_up and hwhm_down will # be different when converting a FWHM from # wavelength to frequency - fwhm = (hwhm_up + hwhm_down) + fwhm = hwhm_up + hwhm_down if quantity == "flux": flux_scaling = wavelength diff --git a/species/read/read_isochrone.py b/species/read/read_isochrone.py index ef2c6015..b1ca24e6 100644 --- a/species/read/read_isochrone.py +++ b/species/read/read_isochrone.py @@ -50,6 +50,7 @@ def __init__( tag: Optional[str] = None, create_regular_grid: bool = False, extrapolate: bool = False, + verbose: bool = True, ) -> None: """ Parameters @@ -76,12 +77,14 @@ def __init__( since there might be inaccuracies in the extrapolated parts of the parameter space, in particular for the cooling curves extracted with - :func:`~species.read.read_isochrone.ReadIsochrone.get_cooling_curve`. + :func:`~species.read.read_isochrone.ReadIsochrone.get_cooling_track`. extrapolate : str DEPRECATED: This parameter has been renamed to ``create_regular_grid`` and will be removed in a future release. Please use the ``create_regular_grid`` parameter instead. + verbose : bool + Print output. Returns ------- @@ -89,14 +92,15 @@ def __init__( None """ - print_section("Read isochrone grid") - self.tag = tag self.extrapolate = extrapolate self.create_regular_grid = create_regular_grid + self.verbose = verbose - print(f"Database tag: {self.tag}") - print(f"Create regular grid: {self.create_regular_grid}") + if self.verbose: + print_section("Read isochrone grid") + print(f"Database tag: {self.tag}") + print(f"Create regular grid: {self.create_regular_grid}") if self.extrapolate: warnings.warn( @@ -191,7 +195,8 @@ def __init__( self.extra_param["feh"] = float(tag_split[1][:-5]) self.extra_param["fsed"] = float(tag_split[2]) - print(f"\nSetting 'extra_param' attribute: {self.extra_param}") + if self.verbose: + print(f"\nSetting 'extra_param' attribute: {self.extra_param}") @typechecked def _read_data( @@ -710,7 +715,7 @@ def get_isochrone( ) @typechecked - def get_cooling_curve( + def get_cooling_track( self, mass: float, ages: Optional[np.ndarray] = None, diff --git a/species/read/read_model.py b/species/read/read_model.py index 27e10d2a..7ab957a6 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -187,11 +187,15 @@ def wavelength_points(self) -> Tuple[np.ndarray, np.ndarray]: wl_index = np.ones(wl_points.shape[0], dtype=bool) else: - wl_index = (wl_points > self.wavel_range[0]) & ( - wl_points < self.wavel_range[1] + wl_index = (wl_points >= self.wavel_range[0]) & ( + wl_points <= self.wavel_range[1] ) index = np.where(wl_index)[0] + # Add extra wavelength points at the boundary to make + # sure that the wavelength range of a filter profile + # is fully included by the model spectrum + if index[0] - 1 >= 0: wl_index[index[0] - 1] = True @@ -233,6 +237,20 @@ def interpolate_grid( grid_points["teff"] <= teff_range[1] ) + # Add extra Teff points at the boundary to make sure + # sure that the Teff prior of a fit is fully included + # in the Teff range that is interpolated + + first_teff = np.where(teff_select)[0][0] + + if first_teff - 1 >= 0: + teff_select[first_teff - 1] = True + + last_teff = np.where(teff_select)[0][-1] + + if last_teff + 1 < teff_select.size: + teff_select[last_teff + 1] = True + grid_points["teff"] = grid_points["teff"][teff_select] else: @@ -261,7 +279,7 @@ def interpolate_grid( grid_points, grid_flux, method="linear", - bounds_error=False, + bounds_error=True, fill_value=np.nan, ) @@ -1836,10 +1854,10 @@ def integrate_spectrum(self, model_param: Dict[str, float]) -> float: applying extinction, the integrated luminosity should in principle be the same as the luminosity calculated directly from the :math:`T_\\mathrm{eff}` and radius parameters. Unless - a particular spectrum from a radiative-convective model had - not fully converged, so it can be useful to check if the - integrated luminosity is indeed consistent with the - :math:`T_\\mathrm{eff}` of the model. + the radiative-convective model had not fully converged for a + particular set of input parameters. It can thus be useful + to check if the integrated luminosity is indeed consistent + with the :math:`T_\\mathrm{eff}` of the model. Parameters ---------- @@ -2071,7 +2089,6 @@ def create_color_color( # with self.open_database() as hdf5_file: # wl_points = np.array(hdf5_file[f"models/{self.model}/wavelength"]) # grid_flux = np.array(hdf5_file[f"models/{self.model}/flux"]) - # print(grid_flux.shape) # # import matplotlib.pyplot as plt # @@ -2100,7 +2117,6 @@ def create_color_color( # # Create list with grid points # # grid_points = list(grid_points.values()) - # print(grid_points) # # # Get the boolean array for selecting the fluxes # # within the requested wavelength range diff --git a/species/util/fit_util.py b/species/util/fit_util.py index 470bfee8..5c602179 100644 --- a/species/util/fit_util.py +++ b/species/util/fit_util.py @@ -296,9 +296,12 @@ def get_residuals( with h5py.File(database_path, "r") as hdf5_file: dset = hdf5_file[f"results/fit/{tag}/samples"] - spectrum = dset.attrs["spectrum"] + if "model_name" in dset.attrs: + model_name = dset.attrs["model_name"] + elif "spectrum" in dset.attrs: + model_name = dset.attrs["spectrum"] binary = dset.attrs["binary"] - print(f"Model: {spectrum}") + print(f"Model: {model_name}") print(f"Binary: {binary}") n_param = dset.attrs["n_param"] @@ -347,7 +350,7 @@ def get_residuals( model_phot = multi_photometry( datatype="model", - spectrum=spectrum, + spectrum=model_name, filters=inc_phot, parameters=parameters, radtrans=radtrans, @@ -378,7 +381,7 @@ def get_residuals( if inc_spec and objectbox.spectrum is not None: res_spec = {} - if spectrum == "petitradtrans": + if model_name == "petitradtrans": # Calculate the petitRADTRANS spectrum only once # Smoothing and resampling not with get_model model = radtrans.get_model(parameters) @@ -393,7 +396,7 @@ def get_residuals( wl_new = objectbox.spectrum[key][0][:, 0] spec_res = objectbox.spectrum[key][3] - if spectrum == "planck": + if model_name == "planck": readmodel = ReadPlanck(wavel_range=wavel_range) model = readmodel.get_spectrum( @@ -411,7 +414,7 @@ def get_residuals( verbose=True, ) - elif spectrum == "petitradtrans": + elif model_name == "petitradtrans": # Smoothing to the instrument resolution flux_smooth = convolve_spectrum( model.wavelength, model.flux, spec_res @@ -431,7 +434,7 @@ def get_residuals( # Resampling to the new wavelength points # is done by the get_model method - readmodel = ReadModel(spectrum, wavel_range=wavel_range) + readmodel = ReadModel(model_name, wavel_range=wavel_range) if "teff_0" in parameters and "teff_1" in parameters: # Binary system @@ -466,7 +469,7 @@ def get_residuals( model_spec = create_box( boxtype="model", - model=spectrum, + model=model_name, wavelength=wl_new, flux=flux_comb, parameters=parameters,