Skip to content

Commit

Permalink
Fixing issue with flux ratio prior, enabling plotting of binary masse…
Browse files Browse the repository at this point in the history 
…s with inc_mass and inc_log_mass
  • Loading branch information
@tomasstolker
tomasstolker committed May 3, 2024
1 parent b211156 commit 01f917c
Show file tree
Hide file tree
Showing 5 changed files with 212 additions and 99 deletions.
73 changes: 37 additions & 36 deletions species/fit/fit_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {}
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
82 changes: 81 additions & 1 deletion species/plot/plot_mcmc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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'."
Expand All @@ -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":
Expand Down
25 changes: 14 additions & 11 deletions species/read/read_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])
)
Expand All @@ -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])
)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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])
)
Expand All @@ -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])
)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 01f917c

Please sign in to comment.