Skip to content

Commit

Permalink
Support for duplicate filter names in FitModel and FitSpectrum
Browse files Browse the repository at this point in the history
  • Loading branch information
Tomas Stolker committed Apr 30, 2020
1 parent 7ef1829 commit 8d20c55
Show file tree
Hide file tree
Showing 6 changed files with 55 additions and 41 deletions.
22 changes: 17 additions & 5 deletions species/analysis/fit_model.py
View file
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
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
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
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
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
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

0 comments on commit 8d20c55

Please sign in to comment.