diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py index 2ec0f97..6f6a75a 100644 --- a/species/fit/fit_model.py +++ b/species/fit/fit_model.py @@ -1049,24 +1049,22 @@ def __init__( # Include prior flux ratio when fitting - self.flux_ratio = [] + self.flux_ratio = {} for param_item in self.bounds: if param_item[:6] == "ratio_": - self.modelpar.append(param_item) print(f"Interpolating {param_item[6:]}...", end="", flush=True) read_model = ReadModel(self.model, filter_name=param_item[6:]) read_model.interpolate_grid(wavel_resample=None, spec_res=None) - self.flux_ratio.append(read_model) + self.flux_ratio[param_item[6:]] = read_model print(" [DONE]") for param_item in self.normal_prior: if param_item[:6] == "ratio_": - self.modelpar.append(param_item) print(f"Interpolating {param_item[6:]}...", end="", flush=True) read_model = ReadModel(self.model, filter_name=param_item[6:]) read_model.interpolate_grid(wavel_resample=None, spec_res=None) - self.flux_ratio.append(read_model) + self.flux_ratio[param_item[6:]] = read_model print(" [DONE]") self.fix_param = {} @@ -1564,6 +1562,40 @@ def _lnlike_func( "so the mass prior can not be applied." ) + elif key[:6] == "ratio_": + filt_name = key[6:] + + param_0 = binary_to_single(param_dict, 0) + phot_flux_0 = self.flux_ratio[filt_name].spectrum_interp( + list(param_0.values()) + )[0][0] + + param_1 = binary_to_single(param_dict, 1) + phot_flux_1 = self.flux_ratio[filt_name].spectrum_interp( + list(param_1.values()) + )[0][0] + + # Uniform prior for the flux ratio + + if f"ratio_{filt_name}" in self.bounds: + ratio_prior = self.bounds[f"ratio_{filt_name}"] + + if ratio_prior[0] > phot_flux_1 / phot_flux_0: + return -np.inf + elif ratio_prior[1] < phot_flux_1 / phot_flux_0: + return -np.inf + + # Normal prior for the flux ratio + + if f"ratio_{filt_name}" in self.normal_prior: + ratio_prior = self.normal_prior[f"ratio_{filt_name}"] + + ln_like += ( + -0.5 + * (phot_flux_1 / phot_flux_0 - ratio_prior[0]) ** 2 + / ratio_prior[1] ** 2 + ) + else: ln_like += ( -0.5 @@ -1589,37 +1621,6 @@ def _lnlike_func( dust_param["powerlaw_ext"] / cross_tmp / 2.5 / np.log10(np.exp(1.0)) ) - # Optionally check the flux ratio of a requested filter - - for filter_idx, model_item in enumerate(self.flux_ratio): - filt_name = model_item.filter_name - - param_0 = binary_to_single(param_dict, 0) - phot_flux_0 = model_item.spectrum_interp(list(param_0.values()))[0][0] - - param_1 = binary_to_single(param_dict, 1) - phot_flux_1 = model_item.spectrum_interp(list(param_1.values()))[0][0] - - # Uniform prior for the flux ratio - - if f"ratio_{filt_name}" in self.bounds: - ratio_prior = self.bounds[f"ratio_{filt_name}"] - if ratio_prior[0] > phot_flux_1 / phot_flux_0: - return -np.inf - elif ratio_prior[1] < phot_flux_1 / phot_flux_0: - return -np.inf - - # Normal prior for the flux ratio - - if f"ratio_{filt_name}" in self.normal_prior: - ratio_prior = self.normal_prior[f"ratio_{filt_name}"] - - ln_like += ( - -0.5 - * (phot_flux_1 / phot_flux_0 - ratio_prior[0]) ** 2 - / ratio_prior[1] ** 2 - ) - for i, obj_item in enumerate(self.objphot): # Get filter name phot_filter = self.modelphot[i].filter_name diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py index b0ee91d..b44f27b 100644 --- a/species/plot/plot_mcmc.py +++ b/species/plot/plot_mcmc.py @@ -680,6 +680,8 @@ def plot_posterior( # Include the derived mass if inc_mass: + check_param = False + if "logg" in box.parameters and "radius" in box.parameters: logg_index = np.argwhere(np.array(box.parameters) == "logg")[0] radius_index = np.argwhere(np.array(box.parameters) == "radius")[0] @@ -693,12 +695,46 @@ def plot_posterior( box.parameters.append("mass") ndim += 1 - else: + check_param = True + + if "logg_0" in box.parameters and "radius_0" in box.parameters: + logg_index = np.argwhere(np.array(box.parameters) == "logg_0")[0] + radius_index = np.argwhere(np.array(box.parameters) == "radius_0")[0] + + mass_samples = logg_to_mass( + samples[..., logg_index], samples[..., radius_index] + ) + + samples = np.append(samples, mass_samples, axis=-1) + + box.parameters.append("mass_0") + ndim += 1 + + check_param = True + + if "logg_1" in box.parameters and "radius_1" in box.parameters: + logg_index = np.argwhere(np.array(box.parameters) == "logg_1")[0] + radius_index = np.argwhere(np.array(box.parameters) == "radius_1")[0] + + mass_samples = logg_to_mass( + samples[..., logg_index], samples[..., radius_index] + ) + + samples = np.append(samples, mass_samples, axis=-1) + + box.parameters.append("mass_1") + ndim += 1 + + check_param = True + + if not check_param: warnings.warn( "Samples with the log(g) and radius are required for 'inc_mass=True'." ) if inc_log_mass: + check_param = False + if "logg" in box.parameters and "radius" in box.parameters: logg_index = np.argwhere(np.array(box.parameters) == "logg")[0] radius_index = np.argwhere(np.array(box.parameters) == "radius")[0] @@ -713,6 +749,40 @@ def plot_posterior( box.parameters.append("log_mass") ndim += 1 + check_param = True + + if "logg_0" in box.parameters and "radius_0" in box.parameters: + logg_index = np.argwhere(np.array(box.parameters) == "logg_0")[0] + radius_index = np.argwhere(np.array(box.parameters) == "radius_0")[0] + + mass_samples = logg_to_mass( + samples[..., logg_index], samples[..., radius_index] + ) + + mass_samples = np.log10(mass_samples) + samples = np.append(samples, mass_samples, axis=-1) + + box.parameters.append("log_mass_0") + ndim += 1 + + check_param = True + + if "logg_1" in box.parameters and "radius_1" in box.parameters: + logg_index = np.argwhere(np.array(box.parameters) == "logg_1")[0] + radius_index = np.argwhere(np.array(box.parameters) == "radius_1")[0] + + mass_samples = logg_to_mass( + samples[..., logg_index], samples[..., radius_index] + ) + + mass_samples = np.log10(mass_samples) + samples = np.append(samples, mass_samples, axis=-1) + + box.parameters.append("log_mass_1") + ndim += 1 + + check_param = True + else: warnings.warn( "Samples with the log(g) and radius are required for 'inc_log_mass=True'." @@ -738,6 +808,16 @@ def plot_posterior( if object_type == "star": samples[:, mass_index] *= constants.M_JUP / constants.M_SUN + if "mass_0" in box.parameters: + mass_index = np.argwhere(np.array(box.parameters) == "mass_0")[0] + if object_type == "star": + samples[:, mass_index] *= constants.M_JUP / constants.M_SUN + + if "mass_1" in box.parameters: + mass_index = np.argwhere(np.array(box.parameters) == "mass_1")[0] + if object_type == "star": + samples[:, mass_index] *= constants.M_JUP / constants.M_SUN + if "log_mass" in box.parameters: mass_index = np.argwhere(np.array(box.parameters) == "log_mass")[0] if object_type == "star": diff --git a/species/read/read_model.py b/species/read/read_model.py index d2e352c..bbc5ba5 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -894,13 +894,15 @@ def get_model( else: for disk_idx in range(100): - if f"disk_teff_{disk_idx}" in model_param and f"disk_radius_{disk_idx}" in model_param: + if ( + f"disk_teff_{disk_idx}" in model_param + and f"disk_radius_{disk_idx}" in model_param + ): n_disk += 1 else: break if n_disk == 1: - readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) @@ -924,7 +926,6 @@ def get_model( flux += flux_interp(self.wl_points) elif n_disk > 1: - readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) @@ -1176,7 +1177,6 @@ def get_model( # Add the blackbody disk components to the luminosity if n_disk == 1: - model_box.parameters["luminosity"] += ( 4.0 * np.pi @@ -1187,12 +1187,15 @@ def get_model( ) # (Lsun) elif n_disk > 1: - for disk_idx in range(n_disk): model_box.parameters["luminosity"] += ( 4.0 * np.pi - * (model_box.parameters[f"disk_radius_{disk_idx}"] * constants.R_JUP) ** 2 + * ( + model_box.parameters[f"disk_radius_{disk_idx}"] + * constants.R_JUP + ) + ** 2 * constants.SIGMA_SB * model_box.parameters[f"disk_teff_{disk_idx}"] ** 4.0 / constants.L_SUN @@ -1400,7 +1403,6 @@ def get_data( break if n_disk == 1: - readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) @@ -1424,7 +1426,6 @@ def get_data( flux += flux_interp(wl_points) elif n_disk > 1: - readplanck = ReadPlanck( (0.9 * self.wavel_range[0], 1.1 * self.wavel_range[-1]) ) @@ -1560,7 +1561,6 @@ def get_data( # Add the blackbody disk components to the luminosity if n_disk == 1: - model_box.parameters["luminosity"] += ( 4.0 * np.pi @@ -1571,12 +1571,15 @@ def get_data( ) # (Lsun) elif n_disk > 1: - for disk_idx in range(n_disk): model_box.parameters["luminosity"] += ( 4.0 * np.pi - * (model_box.parameters[f"disk_radius_{disk_idx}"] * constants.R_JUP) ** 2 + * ( + model_box.parameters[f"disk_radius_{disk_idx}"] + * constants.R_JUP + ) + ** 2 * constants.SIGMA_SB * model_box.parameters[f"disk_teff"] ** 4.0 / constants.L_SUN diff --git a/species/read/read_object.py b/species/read/read_object.py index 36c8b4f..53c06bb 100644 --- a/species/read/read_object.py +++ b/species/read/read_object.py @@ -44,8 +44,8 @@ def __init__(self, object_name: str) -> None: self.database = config["species"]["database"] - with h5py.File(self.database, "r") as h5_file: - if f"objects/{self.object_name}" not in h5_file: + with h5py.File(self.database, "r") as hdf5_file: + if f"objects/{self.object_name}" not in hdf5_file: raise ValueError( f"The object '{self.object_name}' is not " f"present in the database." @@ -73,10 +73,12 @@ def list_filters(self, verbose=True) -> List[str]: if verbose: print(f"Available photometric data for {self.object_name}:") - with h5py.File(self.database, "r") as h5_file: - for tel_item in h5_file[f"objects/{self.object_name}"]: + with h5py.File(self.database, "r") as hdf5_file: + for tel_item in hdf5_file[f"objects/{self.object_name}"]: if tel_item not in ["parallax", "distance", "spectrum"]: - for filt_item in h5_file[f"objects/{self.object_name}/{tel_item}"]: + for filt_item in hdf5_file[ + f"objects/{self.object_name}/{tel_item}" + ]: filter_list.append(f"{tel_item}/{filt_item}") if verbose: @@ -105,25 +107,24 @@ def get_photometry(self, filter_name: str) -> np.ndarray: flux (W m-2 um-1), flux error (W m-2 um-1). """ - with h5py.File(self.database, "r") as h5_file: - if filter_name in h5_file[f"objects/{self.object_name}"]: - obj_phot = np.asarray( - h5_file[f"objects/{self.object_name}/{filter_name}"] - ) - - else: + with h5py.File(self.database, "r") as hdf5_file: + if filter_name not in hdf5_file[f"objects/{self.object_name}"]: raise ValueError( f"There is no photometric data of {self.object_name} " f"available with the {filter_name} filter." ) + obj_phot = np.asarray( + hdf5_file[f"objects/{self.object_name}/{filter_name}"] + ) + return obj_phot @typechecked def get_spectrum(self) -> dict: """ - Function for reading the spectra and covariance matrices of the - object. + Function for reading the spectra and covariance matrices + of the selected object. Returns ------- @@ -131,27 +132,27 @@ def get_spectrum(self) -> dict: Dictionary with spectra and covariance matrices. """ - with h5py.File(self.database, "r") as h5_file: - if f"objects/{self.object_name}/spectrum" in h5_file: + with h5py.File(self.database, "r") as hdf5_file: + if f"objects/{self.object_name}/spectrum" in hdf5_file: spectrum = {} - for item in h5_file[f"objects/{self.object_name}/spectrum"]: - data_group = f"objects/{self.object_name}/spectrum/{item}" + for spec_item in hdf5_file[f"objects/{self.object_name}/spectrum"]: + data_group = f"objects/{self.object_name}/spectrum/{spec_item}" - if f"{data_group}/covariance" not in h5_file: - spectrum[item] = ( - np.asarray(h5_file[f"{data_group}/spectrum"]), + if f"{data_group}/covariance" not in hdf5_file: + spectrum[spec_item] = ( + np.asarray(hdf5_file[f"{data_group}/spectrum"]), None, None, - h5_file[f"{data_group}"].attrs["specres"], + hdf5_file[f"{data_group}"].attrs["specres"], ) else: - spectrum[item] = ( - np.asarray(h5_file[f"{data_group}/spectrum"]), - np.asarray(h5_file[f"{data_group}/covariance"]), - np.asarray(h5_file[f"{data_group}/inv_covariance"]), - h5_file[f"{data_group}"].attrs["specres"], + spectrum[spec_item] = ( + np.asarray(hdf5_file[f"{data_group}/spectrum"]), + np.asarray(hdf5_file[f"{data_group}/covariance"]), + np.asarray(hdf5_file[f"{data_group}/inv_covariance"]), + hdf5_file[f"{data_group}"].attrs["specres"], ) else: @@ -172,16 +173,16 @@ def get_distance(self) -> Tuple[float, float]: Uncertainty (pc). """ - with h5py.File(self.database, "r") as h5_file: - if f"objects/{self.object_name}/parallax" in h5_file: - parallax = np.asarray(h5_file[f"objects/{self.object_name}/parallax"]) + with h5py.File(self.database, "r") as hdf5_file: + if f"objects/{self.object_name}/parallax" in hdf5_file: + parallax = np.asarray(hdf5_file[f"objects/{self.object_name}/parallax"]) calc_dist = 1.0 / (parallax[0] * 1e-3) # (pc) dist_plus = 1.0 / ((parallax[0] - parallax[1]) * 1e-3) - calc_dist dist_minus = calc_dist - 1.0 / ((parallax[0] + parallax[1]) * 1e-3) distance = (calc_dist, (dist_plus + dist_minus) / 2.0) - elif f"objects/{self.object_name}/distance" in h5_file: - distance = np.asarray(h5_file[f"objects/{self.object_name}/distance"]) + elif f"objects/{self.object_name}/distance" in hdf5_file: + distance = np.asarray(hdf5_file[f"objects/{self.object_name}/distance"]) else: raise RuntimeError( @@ -206,10 +207,10 @@ def get_parallax(self) -> Tuple[float, float]: Uncertainty (mas). """ - with h5py.File(self.database, "r") as h5_file: - if f"objects/{self.object_name}/parallax" in h5_file: + with h5py.File(self.database, "r") as hdf5_file: + if f"objects/{self.object_name}/parallax" in hdf5_file: obj_parallax = np.asarray( - h5_file[f"objects/{self.object_name}/parallax"] + hdf5_file[f"objects/{self.object_name}/parallax"] ) else: @@ -245,12 +246,12 @@ def get_absmag( Error on the absolute magnitude. """ - with h5py.File(self.database, "r") as h5_file: + with h5py.File(self.database, "r") as hdf5_file: obj_distance = self.get_distance() - if filter_name in h5_file[f"objects/{self.object_name}"]: + if filter_name in hdf5_file[f"objects/{self.object_name}"]: obj_phot = np.asarray( - h5_file[f"objects/{self.object_name}/{filter_name}"] + hdf5_file[f"objects/{self.object_name}/{filter_name}"] ) else: diff --git a/species/util/plot_util.py b/species/util/plot_util.py index 6abed04..7bbf31e 100644 --- a/species/util/plot_util.py +++ b/species/util/plot_util.py @@ -375,12 +375,26 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]: elif object_type == "star": param[index] = r"$M_\ast$ ($M_\mathrm{\odot}$)" - for i, item in enumerate(ascii_lowercase[1:]): - if f"mass_{i}" in param: - index = param.index(f"mass_{i}") - param[index] = rf"$M_\mathrm{{{item}}}$ ($M_\mathrm{{J}}$)" - else: - break + if "mass_0" in param: + index = param.index("mass_0") + if object_type == "planet": + param[index] = r"$M_\mathrm{1}$ ($M_\mathrm{J}$)" + elif object_type == "star": + param[index] = r"$M_\mathrm{1}$ ($M_\mathrm{\odot}$)" + + if "mass_1" in param: + index = param.index("mass_1") + if object_type == "planet": + param[index] = r"$M_\mathrm{2}$ ($M_\mathrm{J}$)" + elif object_type == "star": + param[index] = r"$M_\mathrm{2}$ ($M_\mathrm{\odot}$)" + + # for i, item in enumerate(ascii_lowercase[1:]): + # if f"mass_{i}" in param: + # index = param.index(f"mass_{i}") + # param[index] = rf"$M_\mathrm{{{item}}}$ ($M_\mathrm{{J}}$)" + # else: + # break if "log_mass" in param: index = param.index("log_mass") @@ -389,17 +403,31 @@ def update_labels(param: List[str], object_type: str = "planet") -> List[str]: elif object_type == "star": param[index] = r"$\log\,M_\ast/M_\mathrm{\odot}$" + if "log_mass_0" in param: + index = param.index("log_mass_0") + if object_type == "planet": + param[index] = r"$\log\,M_\mathrm{1}/M_\mathrm{J}$" + elif object_type == "star": + param[index] = r"$\log\,M_\mathrm{1}/M_\mathrm{\odot}$" + + if "log_mass_1" in param: + index = param.index("log_mass_1") + if object_type == "planet": + param[index] = r"$\log\,M_\mathrm{2}/M_\mathrm{J}$" + elif object_type == "star": + param[index] = r"$\log\,M_\mathrm{2}/M_\mathrm{\odot}$" + if "age" in param: index = param.index("age") param[index] = "Age (Myr)" - if "mass_1" in param: - index = param.index("mass_1") - param[index] = r"$M_\mathrm{b}$ ($M_\mathrm{J}$)" - - if "mass_2" in param: - index = param.index("mass_2") - param[index] = r"$M_\mathrm{c}$ ($M_\mathrm{J}$)" + # if "mass_1" in param: + # index = param.index("mass_1") + # param[index] = r"$M_\mathrm{b}$ ($M_\mathrm{J}$)" + # + # if "mass_2" in param: + # index = param.index("mass_2") + # param[index] = r"$M_\mathrm{c}$ ($M_\mathrm{J}$)" if "entropy" in param: index = param.index("entropy")