Skip to content

Commit

Permalink
Support of FitPlanck for multiple Planck functions, improved document…
Browse files Browse the repository at this point in the history 
…ation
  • Loading branch information
Tomas Stolker committed Jan 20, 2020
1 parent 5d62ec6 commit aef54ce
Show file tree
Hide file tree
Showing 14 changed files with 283 additions and 191 deletions.
2 changes: 1 addition & 1 deletion docs/running.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ Running species

To get a first impression, this page shows what a typical workflow with `species` looks like.

First, the atmospheric model spectra will be downloaded and added to the database. Available photometry of PZ Tel B will be read from :class:`~species.data.companions` dictionary and added to the database. Alternatively, the :func:`~species.data.database.Database.add_object` function can be used for manually creating an planet or other type of object in the database. Then, the model spectra are read from the database, interpolated at the provided parameter values, and stored in a :class:`~species.core.box.ModelBox`. The data of PZ Tel B is read from the database and stored in a :class:`~species.core.box.ObjectBox`. The :class:`~species.core.box.Box` objects are given to the :func:`~species.plot.plot_spectrum.plot_spectrum` function, together with the filter names for which the photometry is available.
First, the atmospheric model spectra will be downloaded and added to the database. Available photometry of PZ Tel B will be read from the :class:`~species.data.companions` dictionary and added to the database. Alternatively, the :func:`~species.data.database.Database.add_object` function can be used for manually creating an planet or other type of object in the database. Then, the model spectra are read from the database, interpolated at the provided parameter values, and stored in a :class:`~species.core.box.ModelBox`. The data of PZ Tel B is read from the database and stored in a :class:`~species.core.box.ObjectBox`. The :class:`~species.core.box.Box` objects are given to the :func:`~species.plot.plot_spectrum.plot_spectrum` function, together with the filter names for which the photometry is available.

The following code can be executed from the command line or within a `Jupyter notebook <https://jupyter.org/>`_.

Expand Down
66 changes: 14 additions & 52 deletions species/analysis/fit_model.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
Module for fitting atmospheric models.
Module with functionalities for fitting atmospheric model specra.
"""

import sys
Expand All @@ -22,7 +22,7 @@ def lnprior(param,
modelpar,
prior=None):
"""
Function for the prior probability.
Internal function for the prior probability.
Parameters
----------
Expand Down Expand Up @@ -76,11 +76,10 @@ def lnlike(param,
synphot,
distance,
spectrum,
instrument,
modelspec,
weighting):
"""
Function for the likelihood probability.
Internal function for the likelihood probability.
Parameters
----------
Expand All @@ -95,8 +94,6 @@ def lnlike(param,
Distance (pc).
spectrum : numpy.ndarray
Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1).
instrument : str
Instrument that was used for the spectrum (currently only 'gpi' possible).
modelspec : species.read.read_model.ReadModel
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
Expand Down Expand Up @@ -151,11 +148,10 @@ def lnprob(param,
distance,
prior,
spectrum,
instrument,
modelspec,
weighting):
"""
Function for the posterior probability.
Internal function for the posterior probability.
Parameters
----------
Expand All @@ -176,8 +172,6 @@ def lnprob(param,
used if set to None.
spectrum : numpy.ndarray
Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1).
instrument : str
Instrument that was used for the spectrum (currently only 'gpi' possible).
modelspec : species.read.read_model.ReadModel
weighting : float, None
Weighting applied to the spectrum when calculating the likelihood function in order
Expand Down Expand Up @@ -205,7 +199,6 @@ def lnprob(param,
synphot,
distance,
spectrum,
instrument,
modelspec,
weighting)

Expand All @@ -217,11 +210,11 @@ def lnprob(param,

class FitModel:
"""
Fit atmospheric model spectra to photometric and spectral data.
Class for fitting atmospheric model spectra to photometric data.
"""

def __init__(self,
objname,
object_name,
filters,
model,
bounds,
Expand All @@ -230,7 +223,7 @@ def __init__(self,
"""
Parameters
----------
objname : str
object_name : str
Object name in the database.
filters : tuple(str, )
Filter IDs for which the photometry is selected. All available photometry of the
Expand All @@ -251,7 +244,7 @@ def __init__(self,
None
"""

self.object = read_object.ReadObject(objname)
self.object = read_object.ReadObject(object_name)
self.distance = self.object.get_distance()

self.model = model
Expand All @@ -260,13 +253,8 @@ def __init__(self,
if not inc_phot and not inc_spec:
raise ValueError('No photometric or spectral data has been selected.')

if self.bounds is not None and 'teff' in self.bounds:
teff_bound = self.bounds['teff']
else:
teff_bound = None

if self.bounds is not None:
readmodel = read_model.ReadModel(self.model, None, teff_bound)
readmodel = read_model.ReadModel(self.model)
bounds_grid = readmodel.get_bounds()

for item in bounds_grid:
Expand All @@ -287,11 +275,11 @@ def __init__(self,

if not filters:
species_db = database.Database()
objectbox = species_db.get_object(objname, None)
filters = objectbox.filter
objectbox = species_db.get_object(object_name, None)
filters = objectbox.filters

for item in filters:
readmodel = read_model.ReadModel(self.model, item, teff_bound)
readmodel = read_model.ReadModel(self.model, filter_name=item)
readmodel.interpolate_model()
self.modelphot.append(readmodel)

Expand All @@ -309,7 +297,7 @@ 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, (0.9, 2.5), teff_bound)
self.modelspec = read_model.ReadModel(self.model, wavel_range=(0.9, 2.5))

else:
self.spectrum = None
Expand Down Expand Up @@ -374,7 +362,7 @@ def run_mcmc(self,
high=self.bounds[item][1],
size=nwalkers)

with Pool(processes=cpu_count()) as pool:
with Pool(processes=cpu_count()):
sampler = emcee.EnsembleSampler(nwalkers,
ndim,
lnprob,
Expand All @@ -386,37 +374,11 @@ def run_mcmc(self,
self.distance,
prior,
self.spectrum,
self.instrument,
self.modelspec,
weighting]))

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

# sampler = emcee.EnsembleSampler(nwalkers=nwalkers,
# dim=ndim,
# lnpostfn=lnprob,
# a=2.,
# args=([self.bounds,
# self.modelpar,
# self.modelphot,
# self.objphot,
# self.synphot,
# self.distance,
# prior,
# self.spectrum,
# self.instrument,
# self.modelspec,
# weighting]))

# progbar = progress.bar.Bar('\rRunning MCMC...',
# max=nsteps,
# suffix='%(percent)d%%')

# for i, _ in enumerate(sampler.sample(initial, iterations=nsteps)):
# progbar.next()

# progbar.finish()

species_db = database.Database()

species_db.add_samples(sampler=sampler,
Expand Down
Loading

0 comments on commit aef54ce

Please sign in to comment.