Skip to content

Commit

Permalink
Improved spectral resampling with spectres, functionalities for inclu…
Browse files Browse the repository at this point in the history 
…ding multiple spectra/instruments and covariance matrices in FitModel
  • Loading branch information
Tomas Stolker committed Feb 3, 2020
1 parent 07be94f commit 87ac0cf
Show file tree
Hide file tree
Showing 21 changed files with 487 additions and 322 deletions.
2 changes: 1 addition & 1 deletion README.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ Documentation can be found at `http://species.readthedocs.io <http://species.rea
Attribution
-----------

Please cite `Stolker et al. (2020) <https://ui.adsabs.harvard.edu/abs/2019arXiv191213316S/>`_ whenever results from *species* are used in a publication. Please also make sure to give credit to the relevant papers regarding the use of the publicly available data that *species* benefits from. The *species* logo is `available <https://people.phys.ethz.ch/~stolkert/species/species_logo.zip>`_ for use in presentations.
Please cite `Stolker et al. (2020) <https://ui.adsabs.harvard.edu/abs/2019arXiv191213316S/abstract/>`_ whenever results from *species* are used in a publication. Please also make sure to give credit to the relevant papers regarding the use of the publicly available data that *species* benefits from. The *species* logo is `available <https://people.phys.ethz.ch/~stolkert/species/species_logo.zip>`_ for use in presentations.

Contributing
------------
Expand Down
2 changes: 1 addition & 1 deletion docs/about.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,4 @@ Questions & feedback
Attribution
-----------

Please cite `Stolker et al. (2020) <https://ui.adsabs.harvard.edu/abs/2019arXiv191213316S/>`_ whenever results from *species* are used in a publication. Please also make sure to give credit the relevant papers regarding the use of the publicly available data that *species* is using, such as the photometric and spectral libraries, atmospheric models, evolutionary models, and photometry of individual objects. The *species* logo is `available <https://people.phys.ethz.ch/~stolkert/species/species_logo.zip>`_ for use in presentations.
Please cite `Stolker et al. (2020) <https://ui.adsabs.harvard.edu/abs/2019arXiv191213316S/abstract>`_ whenever results from *species* are used in a publication. Please also make sure to give credit the relevant papers regarding the use of the publicly available data that *species* is using, such as the photometric and spectral libraries, atmospheric models, evolutionary models, and photometry of individual objects. The *species* logo is `available <https://people.phys.ethz.ch/~stolkert/species/species_logo.zip>`_ for use in presentations.
79 changes: 31 additions & 48 deletions species/analysis/fit_model.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
Module with functionalities for fitting atmospheric model specra.
Module with functionalities for fitting atmospheric model spectra.
"""

import math
Expand Down Expand Up @@ -75,8 +75,7 @@ def lnlike(param,
synphot,
distance,
spectrum,
modelspec,
weighting):
modelspec):
"""
Internal function for the likelihood probability.
Expand All @@ -86,20 +85,14 @@ def lnlike(param,
Parameter values.
modelpar : list(str, )
Parameter names.
modelphot : list('species.read.read_model.ReadModel, )
modelphot : list(species.read.read_model.ReadModel, )
objphot : list(tuple(float, float), )
synphot : list(species.analysis.photometry.SyntheticPhotometry, )
distance : float
Distance (pc).
spectrum : numpy.ndarray
spectrum : dict
Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1).
modelspec : species.read.read_model.ReadModel
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
modelspec : list(species.read.read_model.ReadModel, )
Returns
-------
Expand All @@ -121,19 +114,23 @@ def lnlike(param,
chisq += (item[0]-flux)**2 / item[1]**2

if spectrum is not None:
model = modelspec.get_model(paramdict)
for i, item in enumerate(spectrum.keys()):
model = modelspec[i].get_model(paramdict)

flux_new = spectres.spectres(new_spec_wavs=spectrum[:, 0],
old_spec_wavs=model.wavelength,
spec_fluxes=model.flux,
spec_errs=None)
flux_new = spectres.spectres(new_spec_wavs=spectrum[item][0][:, 0],
old_spec_wavs=model.wavelength,
spec_fluxes=model.flux,
spec_errs=None)

if weighting is None:
chisq += np.nansum((spectrum[:, 1] - flux_new)**2/spectrum[:, 2]**2)
if spectrum[item][1] is not None:
spec_res = spectrum[item][0][:, 1] - flux_new

else:
chisq += (weighting/float(spectrum[:, 0].size)) * \
np.nansum((spectrum[:, 1] - flux_new)**2/spectrum[:, 2]**2)
dot_tmp = np.dot(np.transpose(spec_res), np.linalg.inv(spectrum[item][1]))
chisq += np.dot(dot_tmp, spec_res)

else:
chisq += np.nansum((spectrum[item][0][:, 1] - flux_new)**2 /
spectrum[item][0][:, 2]**2)

return -0.5*chisq

Expand All @@ -147,8 +144,7 @@ def lnprob(param,
distance,
prior,
spectrum,
modelspec,
weighting):
modelspec):
"""
Internal function for the posterior probability.
Expand All @@ -169,15 +165,9 @@ def lnprob(param,
Gaussian prior on one of the parameters. Currently only possible for the mass, e.g.
('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. Not
used if set to None.
spectrum : numpy.ndarray
spectrum : dict
Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1).
modelspec : species.read.read_model.ReadModel
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
modelspec : list(species.read.read_model.ReadModel, )
Returns
-------
Expand All @@ -198,8 +188,7 @@ def lnprob(param,
synphot,
distance,
spectrum,
modelspec,
weighting)
modelspec)

if np.isnan(ln_prob):
ln_prob = -np.inf
Expand Down Expand Up @@ -295,12 +284,14 @@ def __init__(self,

if inc_spec:
self.spectrum = self.object.get_spectrum()
self.instrument = self.object.get_instrument()
self.modelspec = read_model.ReadModel(self.model, wavel_range=(0.9, 2.5))

self.modelspec = []
for key, value in self.spectrum.items():
wavel_range = (0.9*value[0][0, 0], 1.1*value[0][-1, 0])
self.modelspec.append(read_model.ReadModel(self.model, wavel_range=wavel_range))

else:
self.spectrum = None
self.instrument = None
self.modelspec = None

self.modelpar = readmodel.get_parameters()
Expand All @@ -311,8 +302,7 @@ def run_mcmc(self,
nsteps,
guess,
tag,
prior=None,
weighting=None):
prior=None):
"""
Function to run the MCMC sampler.
Expand All @@ -331,20 +321,14 @@ def run_mcmc(self,
Gaussian prior on one of the parameters. Currently only possible for the mass, e.g.
('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. Not
used if set to None.
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
Returns
-------
NoneType
None
"""

sigma = {'teff': 5., 'logg': 0.01, 'feh': 0.01, 'radius': 0.01}
sigma = {'teff': 5., 'logg': 0.01, 'feh': 0.01, 'fsed': 0.01, 'co': 0.01, 'radius': 0.01}

print('Running MCMC...')

Expand Down Expand Up @@ -372,8 +356,7 @@ def run_mcmc(self,
self.distance,
prior,
self.spectrum,
self.modelspec,
weighting]))
self.modelspec]))

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

Expand Down
63 changes: 21 additions & 42 deletions species/analysis/fit_planck.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,7 @@ def lnlike(param,
objphot,
synphot,
distance,
spectrum,
weighting):
spectrum):
"""
Internal function for the likelihood probability.
Expand All @@ -76,12 +75,6 @@ def lnlike(param,
spectrum : numpy.ndarray, None
Spectrum array with the wavelength (micron), flux (W m-2 micron-1), and error
(W m-2 micron-1). Not used if set to None.
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
Returns
-------
Expand All @@ -105,21 +98,26 @@ def lnlike(param,
chisq += (obj_item[0]-flux)**2 / obj_item[1]**2

if spectrum is not None:
readplanck = read_planck.ReadPlanck((0.9*spectrum[0, 0], 1.1*spectrum[-1, 0]))
for i, item in enumerate(spectrum.keys()):
readplanck = read_planck.ReadPlanck((0.9*spectrum[item][0][0, 0],
1.1*spectrum[item][0][-1, 0]))

model = readplanck.get_spectrum(paramdict, 100.)
model = readplanck.get_spectrum(paramdict, 100.)

flux_new = spectres.spectres(new_spec_wavs=spectrum[:, 0],
old_spec_wavs=model.wavelength,
spec_fluxes=model.flux,
spec_errs=None)
flux_new = spectres.spectres(new_spec_wavs=spectrum[item][0][:, 0],
old_spec_wavs=model.wavelength,
spec_fluxes=model.flux,
spec_errs=None)

if weighting is None:
chisq += np.nansum((spectrum[:, 1] - flux_new)**2/spectrum[:, 2]**2)
if spectrum[item][1] is not None:
spec_res = spectrum[item][0][:, 1] - flux_new

else:
chisq += (weighting/float(spectrum[:, 0].size)) * \
np.nansum((spectrum[:, 1] - flux_new)**2/spectrum[:, 2]**2)
dot_tmp = np.dot(np.transpose(spec_res), np.linalg.inv(spectrum[item][1]))
chisq += np.dot(dot_tmp, spec_res)

else:
chisq += np.nansum((spectrum[item][0][:, 1] - flux_new)**2 /
spectrum[item][0][:, 2]**2)

return -0.5*chisq

Expand All @@ -129,8 +127,7 @@ def lnprob(param,
objphot,
synphot,
distance,
spectrum,
weighting):
spectrum):
"""
Internal function for the posterior probability.
Expand All @@ -152,12 +149,6 @@ def lnprob(param,
spectrum : numpy.ndarray, None
Spectrum array with the wavelength (micron), flux (W m-2 micron-1), and error
(W m-2 micron-1). Not used if set to None.
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
Returns
-------
Expand All @@ -176,8 +167,7 @@ def lnprob(param,
objphot,
synphot,
distance,
spectrum,
weighting)
spectrum)

if np.isnan(ln_prob):
ln_prob = -np.inf
Expand Down Expand Up @@ -267,18 +257,14 @@ def __init__(self,

if inc_spec:
self.spectrum = self.object.get_spectrum()
self.instrument = self.object.get_instrument()

else:
self.spectrum = None
self.instrument = None

def run_mcmc(self,
nwalkers,
nsteps,
guess,
tag,
weighting=None):
tag):
"""
Function to run the MCMC sampler.
Expand All @@ -295,12 +281,6 @@ def run_mcmc(self,
``{'teff': [1500., 1000.], 'radius': [1., 2.]``.
tag : str
Database tag where the MCMC samples are stored.
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
to not have a spectrum dominate the chi-squared value. For example, with `weighting=3`
then all combined spectrum points (e.g. covering the YJH bandpasses) have a weighted
that is equal to three photometry points. The spectrum data points have an equal
weighting as the photometry points if set to None.
Returns
-------
Expand Down Expand Up @@ -348,8 +328,7 @@ def run_mcmc(self,
self.objphot,
self.synphot,
self.distance,
self.spectrum,
weighting]))
self.spectrum]))

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

Expand Down
32 changes: 11 additions & 21 deletions species/data/ames_cond.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,9 @@
import tarfile
import urllib.request

import spectres
import numpy as np

from scipy.interpolate import interp1d

from species.util import data_util


Expand Down Expand Up @@ -84,9 +83,6 @@ def add_ames_cond(input_path,
for filename in sorted(file_list):

if filename.startswith('lte') and filename.endswith('.gz'):
print_message = f'Adding AMES-Cond model spectra... {filename}'
print(f'\r{print_message:<80}', end='')

teff_val = float(filename[3:5])*100.
logg_val = float(filename[6:9])
feh_val = float(filename[10:13])
Expand All @@ -98,6 +94,9 @@ def add_ames_cond(input_path,
if teff_val < teff_range[0] or teff_val > teff_range[1]:
continue

print_message = f'Adding AMES-Cond model spectra... {filename}'
print(f'\r{print_message:<80}', end='')

data_wavel = []
data_flux = []

Expand Down Expand Up @@ -164,22 +163,13 @@ def add_ames_cond(input_path,
if np.all(np.diff(data[:, 0]) < 0):
raise ValueError('The wavelengths are not all sorted by increasing value.')

indices = np.where((data[:, 0] >= wavel_range[0]) &
(data[:, 0] <= wavel_range[1]))[0]

if indices.size > 0:
teff.append(teff_val)
logg.append(logg_val)

data = data[indices, :]

flux_interp = interp1d(data[:, 0],
data[:, 1],
kind='linear',
bounds_error=False,
fill_value=1e-100)
teff.append(teff_val)
logg.append(logg_val)

flux.append(flux_interp(wavelength))
try:
flux.append(spectres.spectres(wavelength, data[:, 0], data[:, 1]))
except ValueError:
flux.append(np.zeros(wavelength.shape[0]))

data_sorted = data_util.sort_data(np.asarray(teff),
np.asarray(logg),
Expand All @@ -189,7 +179,7 @@ def add_ames_cond(input_path,
wavelength,
np.asarray(flux))

data_util.write_data('ames-cond', ('teff', 'logg'), database, data_sorted)
data_util.write_data('ames-cond', ['teff', 'logg'], database, data_sorted)

print_message = 'Adding AMES-Cond model spectra... [DONE]'
print(f'\r{print_message:<80}')
Loading

0 comments on commit 87ac0cf

Please sign in to comment.