Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support for duplicate filter names in FitModel and FitSpectrum #9

Merged
merged 1 commit into from
Apr 30, 2020
Merged

Support for duplicate filter names in FitModel and FitSpectrum #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 17 additions & 5 deletions species/analysis/fit_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,9 +128,15 @@ def lnlike(param,
chisq = 0.

if objphot is not None:
for i, item in enumerate(objphot):
for i, obj_item in enumerate(objphot):
flux = scaling * modelphot[i].spectrum_interp(list(paramdict.values()))
chisq += (item[0]-flux)**2 / item[1]**2

if obj_item.ndim == 1:
chisq += (obj_item[0] - flux)**2 / obj_item[1]**2

else:
for i in range(obj_item.shape[1]):
chisq += (obj_item[0, i] - flux)**2 / obj_item[1, i]**2

if spectrum is not None:
for i, item in enumerate(spectrum.keys()):
Expand Down Expand Up @@ -309,7 +315,7 @@ def __init__(self,
print(f' [DONE]')

obj_phot = self.object.get_photometry(item)
self.objphot.append((obj_phot[2], obj_phot[3]))
self.objphot.append(np.array([obj_phot[2], obj_phot[3]]))

else:
self.objphot = None
Expand Down Expand Up @@ -572,9 +578,15 @@ def lnlike_multinest(cube, n_dim, n_param):
chisq = 0.

if self.objphot is not None:
for i, item in enumerate(self.objphot):
for i, obj_item in enumerate(self.objphot):
flux = scaling * self.modelphot[i].spectrum_interp(list(paramdict.values()))
chisq += (item[0]-flux)**2 / item[1]**2

if obj_item.ndim == 1:
chisq += (obj_item[0] - flux)**2 / obj_item[1]**2

else:
for i in range(obj_item.shape[1]):
chisq += (obj_item[0, i] - flux)**2 / obj_item[1, i]**2

if self.spectrum is not None:
for i, item in enumerate(self.spectrum.keys()):
Expand Down
9 changes: 7 additions & 2 deletions species/analysis/fit_planck.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,12 @@ def lnlike(param,
readplanck = read_planck.ReadPlanck(filter_name=synphot[i].filter_name)
flux = readplanck.get_flux(paramdict, synphot=synphot[i])[0]

chisq += (obj_item[0]-flux)**2 / obj_item[1]**2
if obj_item.ndim == 1:
chisq += (obj_item[0]-flux)**2 / obj_item[1]**2

else:
for i in range(obj_item.shape[1]):
chisq += (obj_item[0, i]-flux)**2 / obj_item[1, i]**2

if spectrum is not None:
for i, item in enumerate(spectrum.keys()):
Expand Down Expand Up @@ -249,7 +254,7 @@ def __init__(self,
self.synphot.append(sphot)

obj_phot = self.object.get_photometry(item)
self.objphot.append((obj_phot[2], obj_phot[3]))
self.objphot.append(np.array([obj_phot[2], obj_phot[3]]))

else:
self.synphot = None
Expand Down
35 changes: 10 additions & 25 deletions species/analysis/fit_spectrum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,16 +18,14 @@ def lnprob(param,
bounds,
modelpar,
objphot,
specphot,
bands):
specphot):
"""
Internal function for the posterior probability.

Parameters
----------
param : numpy.ndarray
Values of the main scaling parameter and optionally additional band-dependent scaling
parameters.
Value of the scaling parameter.
bounds : dict
Boundaries of the main scaling parameter.
modelpar : list(str, )
Expand All @@ -37,9 +35,6 @@ def lnprob(param,
specphot : list(float, )
Synthetic photometry of the calibration spectrum for the same filters as the photometry
of the object.
bands : bool
Use band-dependent scaling parameters in addition to the main scaling parameter which
is used for the full spectrum.

Returns
-------
Expand All @@ -61,12 +56,13 @@ def lnprob(param,

else:
chisq = 0.
for i, _ in enumerate(objphot):
if bands:
chisq += (objphot[i][0] - param[0]*param[i+1]*specphot[i])**2 / objphot[i][1]**2
for i, obj_item in enumerate(objphot):
if obj_item.ndim == 1:
chisq += (obj_item[0] - param[0]*specphot[i])**2 / obj_item[1]**2

else:
chisq += (objphot[i][0] - param[0]*specphot[i])**2 / objphot[i][1]**2
for i in range(obj_item.shape[1]):
chisq += (obj_item[0, i] - param[0]*specphot[i])**2 / obj_item[1, i]**2

ln_prob = ln_prior - 0.5*chisq

Expand Down Expand Up @@ -128,16 +124,15 @@ def __init__(self,
self.specphot.append(spec_phot[0])

obj_phot = self.object.get_photometry(item)
self.objphot.append((obj_phot[2], obj_phot[3]))
self.objphot.append(np.array([obj_phot[2], obj_phot[3]]))

self.modelpar = ['scaling']

def run_mcmc(self,
nwalkers,
nsteps,
guess,
tag,
bands=False):
tag):
"""
Function to run the MCMC sampler.

Expand All @@ -151,9 +146,6 @@ def run_mcmc(self,
Guess of the scaling parameter.
tag : str
Database tag where the MCMC samples are stored.
bands : bool
Use band-dependent scaling parameters in addition to the main scaling parameter which
is used for the full spectrum.

Returns
-------
Expand All @@ -162,12 +154,6 @@ def run_mcmc(self,

print('Running MCMC...')

if bands:
ndim = 1 + len(self.objphot)

else:
ndim = 1

initial = np.zeros((nwalkers, ndim))
initial[:, 0] = guess['scaling'] + np.random.normal(0, 1e-1*guess['scaling'], nwalkers)

Expand All @@ -184,8 +170,7 @@ def run_mcmc(self,
args=([self.bounds,
self.modelpar,
self.objphot,
self.specphot,
bands]))
self.specphot]))

ens_sampler.run_mcmc(initial, nsteps, progress=True)

Expand Down
10 changes: 6 additions & 4 deletions species/data/companions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,14 +102,16 @@ def get_data():
'PDS 70 b': {'distance': (113.43, 0.52),
'app_mag': {'Paranal/SPHERE.IRDIS_D_H23_2': (17.94, 0.24), # Keppler et al. 2018
'Paranal/SPHERE.IRDIS_D_H23_3': (17.95, 0.17), # Keppler et al. 2018
'Paranal/SPHERE.IRDIS_D_K12_1': [(16.78, 0.31), (16.68, 0.04)], # Stolker et al. in prep.
'Paranal/SPHERE.IRDIS_D_K12_2': [(16.23, 0.32), (16.35, 0.07)], # Stolker et al. in prep.
'Paranal/SPHERE.IRDIS_D_K12_1': [(16.78, 0.31), (16.72, 0.05)], # Stolker et al. in prep.
'Paranal/SPHERE.IRDIS_D_K12_2': [(16.23, 0.32), (16.38, 0.06)], # Stolker et al. in prep.
'Paranal/NACO.Lp': (14.08, 0.33), # Stolker et al. in prep.
'Paranal/NACO.NB405': (13.91, 0.34), # Stolker et al. in prep.
'Paranal/NACO.Mp': (13.64, 0.22)}}, # Stolker et al. in prep.
'Paranal/NACO.Mp': (13.64, 0.22), # Stolker et al. in prep.
'Keck/NIRC2.Lp': (14.64, 0.18)}}, # Wang et al. 2020

'PDS 70 c': {'distance': (113.43, 0.52),
'app_mag': {'Paranal/NACO.NB405': (15.05, 0.59)}}, # Stolker et al. in prep.
'app_mag': {'Paranal/NACO.NB405': (15.05, 0.59), # Stolker et al. in prep.
'Keck/NIRC2.Lp': (15.5, 0.46)}}, # Wang et al. 2020

'2M1207 b': {'distance': (64.42, 0.65),
'app_mag': {'HST/NICMOS1.F090M': (22.58, 0.35), # Song et al. 2006
Expand Down
19 changes: 14 additions & 5 deletions species/plot/plot_spectrum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -348,7 +348,7 @@ def plot_spectrum(boxes,
if plot_kwargs[j]:
ax1.plot(wavelength, masked/scaling, zorder=1, **plot_kwargs[j])
else:
ax1.plot(wavelength, masked/scaling, lw=0.2, alpha=0.5, zorder=1)
ax1.plot(wavelength, masked/scaling, color='gray', lw=0.2, alpha=0.5, zorder=1)

elif isinstance(boxitem, box.PhotometryBox):
for i, item in enumerate(boxitem.wavelength):
Expand Down Expand Up @@ -409,8 +409,12 @@ def plot_spectrum(boxes,
mask=np.isnan(boxitem.spectrum[key][0]))

if not plot_kwargs[j] or key not in plot_kwargs[j]:
ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling,
ms=2, marker='s', zorder=2.5, ls='none')
plot_obj = ax1.errorbar(masked[:, 0], masked[:, 1]/scaling,
yerr=masked[:, 2]/scaling, ms=2, marker='s',
zorder=2.5, ls='none')

plot_kwargs[j][key] = {'marker': 's', 'ms': 2., 'ls': 'none',
'color': plot_obj[0].get_color()}

else:
ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling,
Expand Down Expand Up @@ -470,8 +474,13 @@ def plot_spectrum(boxes,

elif residuals.photometry[item].ndim == 2:
for i in range(residuals.photometry[item].shape[1]):
ax3.plot(residuals.photometry[item][0, i], residuals.photometry[item][1, i], zorder=2,
**plot_kwargs[obj_index][item][i])
if isinstance(plot_kwargs[obj_index][item], list):
ax3.plot(residuals.photometry[item][0, i], residuals.photometry[item][1, i], zorder=2,
**plot_kwargs[obj_index][item][i])

else:
ax3.plot(residuals.photometry[item][0, i], residuals.photometry[item][1, i], zorder=2,
**plot_kwargs[obj_index][item])

res_max = np.nanmax(np.abs(residuals.photometry[item][1]))

Expand Down
1 change: 1 addition & 0 deletions species/read/read_object.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ def get_photometry(self,
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:
raise ValueError(f'There is no photometric data of {self.object_name} '
f'available with the {filter_name} filter.')
Expand Down