Skip to content

Commit

Permalink
Added weights parameter in FitModel, included disk luminosity in plot…
Browse files Browse the repository at this point in the history 
…_posterior
  • Loading branch information
@tomasstolker
tomasstolker committed Feb 17, 2021
1 parent f51ea4e commit f4d3221
Show file tree
Hide file tree
Showing 6 changed files with 106 additions and 22 deletions.
63 changes: 49 additions & 14 deletions species/analysis/fit_model.py
View file
Expand Up @@ -383,7 +383,8 @@ def __init__(self,
List[Tuple[float, float]]]],
inc_phot: Union[bool, List[str]] = True,
inc_spec: Union[bool, List[str]] = True,
fit_corr: Optional[List[str]] = None) -> None:
fit_corr: Optional[List[str]] = None,
weights: Optional[Dict[str, float]] = None) -> None:
"""
The grid of spectra is linearly interpolated for each photometric point and spectrum while
taking into account the filter profile, spectral resolution, and wavelength sampling.
Expand Down Expand Up @@ -563,6 +564,12 @@ def __init__(self,
fit_corr : list(str), None
List with spectrum names for which the correlation length and fractional amplitude are
fitted (see Wang et al. 2020).
weights : dict(str, float), None
Weights to be applied to the log-likelihood components of the different spectroscopic
and photometric data that are provided with ``inc_spec`` and ``inc_phot``. This
parameter can for example be used to bias the weighting of the photometric data points.
An equal weighting is applied if the argument is set to ``None``. Only supported by
``run_ultranest`` and ``run_multinest``.
Returns
-------
Expand Down Expand Up @@ -912,6 +919,27 @@ def __init__(self,
for i, item in enumerate(self.modelpar):
self.cube_index[item] = i

# Weighting of the photometric and spectroscopic data

print('Weights for the log-likelihood function:')

if weights is None:
self.weights = {}
else:
self.weights = weights

for item in inc_spec:
if item not in self.weights:
self.weights[item] = 1.

print(f' - {item} = {self.weights[item]:.2e}')

for item in inc_phot:
if item not in self.weights:
self.weights[item] = 1.

print(f' - {item} = {self.weights[item]:.2e}')

@typechecked
def run_mcmc(self,
tag: str,
Expand Down Expand Up @@ -1249,10 +1277,11 @@ def lnlike_func(self,
n_grains = dust_param['powerlaw_ext'] / cross_tmp / 2.5 / np.log10(np.exp(1.))

for i, obj_item in enumerate(self.objphot):
if self.modelphot[i] is not None:
phot_filter = self.modelphot[i].filter_name
else:
phot_filter = None
# Get filter name
phot_filter = self.modelphot[i].filter_name

# Shortcut for weight
weight = self.weights[phot_filter]

if self.model == 'planck':
readplanck = read_planck.ReadPlanck(filter_name=phot_filter)
Expand Down Expand Up @@ -1305,10 +1334,10 @@ def lnlike_func(self,
if self.model == 'powerlaw' and f'{phot_filter}_error' in param_dict:
phot_var += param_dict[f'{phot_filter}_error']**2 * obj_item[0]**2

ln_like += -0.5 * (obj_item[0] - phot_flux)**2 / phot_var
ln_like += -0.5 * weight * (obj_item[0] - phot_flux)**2 / phot_var

# Only required when fitting an error inflation
ln_like += -0.5 * np.log(2.*np.pi*phot_var)
ln_like += -0.5 * weight * np.log(2.*np.pi*phot_var)

else:
phot_var = obj_item[1, j]**2
Expand All @@ -1317,14 +1346,17 @@ def lnlike_func(self,
phot_var += param_dict[f'{phot_filter}_error']**2 * obj_item[0, j]**2

for j in range(obj_item.shape[1]):
ln_like += -0.5 * (obj_item[0, j] - phot_flux)**2 / phot_var
ln_like += -0.5 * weight * (obj_item[0, j] - phot_flux)**2 / phot_var

# Only required when fitting an error inflation
ln_like += -0.5 * np.log(2.*np.pi*phot_var)
ln_like += -0.5 * weight * np.log(2.*np.pi*phot_var)

for i, item in enumerate(self.spectrum.keys()):
# Calculate or interpolate the model spectrum

# Shortcut for the weight
weight = self.weights[item]

if self.model == 'planck':
# Calculate a blackbody spectrum
readplanck = read_planck.ReadPlanck((0.9*self.spectrum[item][0][0, 0],
Expand Down Expand Up @@ -1414,10 +1446,10 @@ def lnlike_func(self,

if self.spectrum[item][2] is not None:
# Use the inverted covariance matrix
ln_like += -0.5 * np.dot(data_flux-model_flux,
ln_like += -0.5 * weight * np.dot(data_flux-model_flux,
np.dot(data_cov_inv, data_flux-model_flux))

ln_like += -0.5 * np.nansum(np.log(2.*np.pi*data_var))
ln_like += -0.5 * weight * np.nansum(np.log(2.*np.pi*data_var))

else:
if item in self.fit_corr:
Expand All @@ -1435,12 +1467,15 @@ def lnlike_func(self,
dot_tmp = np.dot(data_flux-model_flux,
np.dot(np.linalg.inv(cov_matrix), data_flux-model_flux))

ln_like += -0.5*dot_tmp - 0.5*np.nansum(np.log(2.*np.pi*data_var))
ln_like += -0.5 * weight * dot_tmp
ln_like += -0.5 * np.nansum(np.log(2.*np.pi*data_var))

else:
# Calculate the chi-square without a covariance matrix
ln_like += np.nansum(-0.5 * (data_flux-model_flux)**2 / data_var -
0.5 * np.log(2.*np.pi*data_var))
chi_sq = -0.5 * weight * (data_flux-model_flux)**2 / data_var
chi_sq += -0.5 * weight * np.log(2.*np.pi*data_var)

ln_like += np.nansum(chi_sq)

return ln_like

Expand Down
30 changes: 28 additions & 2 deletions species/data/companions.py
View file
Expand Up @@ -2,7 +2,7 @@
Module for extracting data of directly imaged planets and brown dwarfs.
"""

from typing import Union, Dict, Tuple, List
from typing import Dict, List, Tuple, Union

from typeguard import typechecked

Expand Down Expand Up @@ -255,6 +255,32 @@ def get_data() -> Dict[str, Dict[str, Union[Tuple[float, float],
'Paranal/SPHERE.IRDIS_D_H23_3': (13.94, 0.17), # Cheetham et al. 2018
'Paranal/SPHERE.IRDIS_D_K12_1': (13.77, 0.17), # Cheetham et al. 2018
'Paranal/SPHERE.IRDIS_D_K12_2': (13.45, 0.19), # Cheetham et al. 2018
'Paranal/NACO.Lp': (13.09, 0.17)}}} # Cheetham et al. 2018
'Paranal/NACO.Lp': (13.09, 0.17)}}, # Cheetham et al. 2018

'TYC 8988 b': {'distance': (94.6, 0.3),
'app_mag': {'Paranal/SPHERE.IRDIS_D_Y23_2': (17.03, 0.21), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_Y23_3': (16.67, 0.16), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_J23_2': (16.27, 0.08), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_J23_3': (15.73, 0.07), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_H23_2': (15.11, 0.08), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_H23_3': (14.78, 0.07), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_K12_1': (14.44, 0.04), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_D_K12_2': (14.07, 0.04), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_B_J': (15.73, 0.38), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_B_H': (15.87, 0.38), # Bohn et al. 2019
'Paranal/SPHERE.IRDIS_B_Ks': (14.70, 0.14), # Bohn et al. 2019
'Paranal/NACO.Lp': (13.30, 0.08), # Bohn et al. 2019
'Paranal/NACO.Mp': (13.08, 0.20)}}, # Bohn et al. 2019

'TYC 8988 c': {'distance': (94.6, 0.3),
'app_mag': {'Paranal/SPHERE.IRDIS_D_Y23_3': (22.37, 0.31), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_J23_2': (21.81, 0.22), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_J23_3': (21.17, 0.15), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_H23_2': (19.78, 0.08), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_H23_3': (19.32, 0.06), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_K12_1': (18.34, 0.04), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_D_K12_2': (17.85, 0.09), # Bohn et al. 2020
'Paranal/SPHERE.IRDIS_B_H': (19.69, 0.23), # Bohn et al. 2020
'Paranal/NACO.Lp': (16.29, 0.21)}}} # Bohn et al. 2020

return data
19 changes: 14 additions & 5 deletions species/plot/plot_mcmc.py
View file
Expand Up @@ -199,19 +199,28 @@ def plot_posterior(tag: str,
teff_index = np.argwhere(np.array(box.parameters) == 'teff')[0]
radius_index = np.argwhere(np.array(box.parameters) == 'radius')[0]

luminosity = 4. * np.pi * (samples[..., radius_index]*constants.R_JUP)**2 * \
lum_planet = 4. * np.pi * (samples[..., radius_index]*constants.R_JUP)**2 * \
constants.SIGMA_SB * samples[..., teff_index]**4. / constants.L_SUN

if 'disk_teff' in box.parameters and 'disk_radius' in box.parameters:
teff_index = np.argwhere(np.array(box.parameters) == 'disk_teff')[0]
radius_index = np.argwhere(np.array(box.parameters) == 'disk_radius')[0]

luminosity += 4. * np.pi * (samples[..., radius_index]*constants.R_JUP)**2 * \
lum_disk = 4. * np.pi * (samples[..., radius_index]*constants.R_JUP)**2 * \
constants.SIGMA_SB * samples[..., teff_index]**4. / constants.L_SUN

samples = np.append(samples, np.log10(luminosity), axis=-1)
box.parameters.append('luminosity')
ndim += 1
samples = np.append(samples, np.log10(lum_planet+lum_disk), axis=-1)
box.parameters.append('luminosity')
ndim += 1

samples = np.append(samples, lum_disk/lum_planet, axis=-1)
box.parameters.append('luminosity_disk_planet')
ndim += 1

else:
samples = np.append(samples, np.log10(lum_planet), axis=-1)
box.parameters.append('luminosity')
ndim += 1

elif 'teff_0' in box.parameters and 'radius_0' in box.parameters:
luminosity = 0.
Expand Down
10 changes: 10 additions & 0 deletions species/read/read_model.py
View file
Expand Up @@ -797,6 +797,11 @@ def get_model(self,
model_box.parameters['radius'] * constants.R_JUP)**2 * constants.SIGMA_SB * \
model_box.parameters['teff']**4. / constants.L_SUN # (Lsun)

if 'disk_teff' in model_box.parameters and 'disk_radius' in model_box.parameters:
model_box.parameters['luminosity'] += 4. * np.pi * (
model_box.parameters['disk_radius'] * constants.R_JUP)**2 * constants.SIGMA_SB * \
model_box.parameters['disk_teff']**4. / constants.L_SUN # (Lsun)

return model_box

@typechecked
Expand Down Expand Up @@ -926,6 +931,11 @@ def get_data(self,
model_box.parameters['radius'] * constants.R_JUP)**2 * constants.SIGMA_SB * \
model_box.parameters['teff']**4. / constants.L_SUN # (Lsun)

if 'disk_teff' in model_box.parameters and 'disk_radius' in model_box.parameters:
model_box.parameters['luminosity'] += 4. * np.pi * (
model_box.parameters['disk_radius'] * constants.R_JUP)**2 * constants.SIGMA_SB * \
model_box.parameters['disk_teff']**4. / constants.L_SUN # (Lsun)

return model_box

@typechecked
Expand Down
2 changes: 1 addition & 1 deletion species/util/data_util.py
View file
Expand Up @@ -33,7 +33,7 @@ def update_sptype(sptypes: np.ndarray) -> List[str]:

sptypes_updated = []

for i, spt_item in enumerate(sptypes):
for spt_item in sptypes:

if spt_item == 'None':
sptypes_updated.append('None')
Expand Down
4 changes: 4 additions & 0 deletions species/util/plot_util.py
View file
Expand Up @@ -170,6 +170,10 @@ def update_labels(param: List[str]) -> List[str]:
index = param.index('luminosity_ratio')
param[index] = r'$\mathregular{log}\,\mathregular{L_1}/\mathregular{L_2}$'

if 'luminosity_disk_planet' in param:
index = param.index('luminosity_disk_planet')
param[index] = r'$\mathregular{L_{disk}}/\mathregular{L_{atm}}$'

if 'lognorm_radius' in param:
index = param.index('lognorm_radius')
param[index] = r'$\mathregular{log}\,\mathregular{r_g}$'
Expand Down

0 comments on commit f4d3221

Please sign in to comment.