Collaborate with `species`

[kelle]

The species package, maintained by @tomasstolker, does a bunch of analysis stuff on brown dwarf data and also uses a HDF5 database that includes observed spectra. It would be awesome to make SIMPLE and species play nice together!

[tomasstolker]

Hi @kelle!

I have added the add_simple_object function to species, which retrieves SIMPLE data of a single object and stores this in the species database.

The resolving power of a given instrument/mode is required for fitting model spectra. That information is currently not available in SIMPLE I think? Perhaps that is something to consider to include in Instruments.json. Currently, there is a manual input requested for the resolution when using add_simple_object and spectra are available.

Another interesting use case would be retrieving a subset of objects from SIMPLE, e.g. based on spectral type or gravity, to use this for a color-magnitude diagram or empirical spectrum comparison. I will leave that for another time.

Let me know if you have any feedback and/or suggestions!

[tomasstolker]

Hereby an example for fitting a grid of spectra to the photometry and spectra:

from species import SpeciesInit
from species.data.database import Database
from species.fit.fit_model import FitModel
from species.read.read_model import ReadModel
from species.plot.plot_mcmc import plot_posterior
from species.plot.plot_spectrum import plot_spectrum
from species.util.box_util import update_objectbox
from species.util.fit_util import get_residuals, multi_photometry

object_name = 'Trappist-1'
tag = 'testing'
model_name = 'sphinx'

inc_phot = ['2MASS/2MASS.J', '2MASS/2MASS.H', '2MASS/2MASS.Ks']
inc_spec = ['IRTF SpeX Prism']

SpeciesInit()

database = Database()

database.add_model('sphinx', teff_range=(2000., 3000.))

database.add_simple_object(object_name, simple_database=None)

fit = FitModel(object_name=object_name,
               model=model_name,
               bounds={'radius': (1., 10.),
                       '2MASS/2MASS_error': (0., 20.),
                       'IRTF SpeX Prism': (None, (0., 10.), None)},
               inc_phot=inc_phot,
               inc_spec=inc_spec,
               fit_corr=None,
               normal_prior=None)

fit.run_dynesty(tag=tag,
                n_live_points=200,
                resume=False,
                output='./dynesty/',
                evidence_tolerance=1.0)

plot_posterior(tag=tag,
               burnin=None,
               title=None,
               offset=(-0.3, -0.3),
               inc_luminosity=True,
               inc_mass=False,
               object_type='star',
               param_inc=None,
               output='posterior.png')

samples = database.get_mcmc_spectra(tag=tag,
                                    burnin=None,
                                    random=50,
                                    wavel_range=(0.2, 20.),
                                    spec_res=100.)

best = database.get_median_sample(tag)
print(best)

read_model = ReadModel(model=model_name, wavel_range=(0.2, 20.))
model_box = read_model.get_model(model_param=best, spec_res=100.)

object_box = database.get_object(object_name=object_name,
                                 inc_phot=inc_phot,
                                 inc_spec=inc_spec)

object_box = update_objectbox(object_box, best, model=model_name)

residuals = get_residuals(datatype='model',
                          spectrum=model_name,
                          parameters=best,
                          objectbox=object_box,
                          inc_phot=inc_phot,
                          inc_spec=inc_spec)

synphot = multi_photometry(datatype='model',
                           spectrum=model_name,
                           filters=object_box.filters,
                           parameters=best)

plot_spectrum(boxes=[samples, model_box, object_box],
              residuals=residuals,
              plot_kwargs=[{'ls': '-', 'lw': 0.2, 'color': 'gray'},
                           {'ls': '-', 'lw': 1., 'color': 'black'},
                           {'GAIA/GAIA3.G': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:green', 'mec': 'black', 'ls': 'none', 'label': 'GAIA'},
                            'GAIA/GAIA3.Gbp': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:green', 'mec': 'black', 'ls': 'none'},
                            'GAIA/GAIA3.Grp': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:green', 'mec': 'black', 'ls': 'none'},
                            'GAIA/GAIA3.Grvs': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:green', 'mec': 'black', 'ls': 'none'},
                            '2MASS/2MASS.J': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:orange', 'mec': 'black', 'ls': 'none', 'label': '2MASS'},
                            '2MASS/2MASS.H': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:orange', 'mec': 'black', 'ls': 'none'},
                            '2MASS/2MASS.Ks': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:orange', 'mec': 'black', 'ls': 'none'},
                            'WISE/WISE.W1': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:blue', 'mec': 'black', 'ls': 'none', 'label': 'WISE'},
                            'WISE/WISE.W2': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:blue', 'mec': 'black', 'ls': 'none'},
                            'WISE/WISE.W3': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:blue', 'mec': 'black', 'ls': 'none'},
                            'WISE/WISE.W4': {'marker': 's', 'ms': 4, 'mew': 1., 'color': 'tab:blue', 'mec': 'black', 'ls': 'none'},
                            'IRTF SpeX Prism': {'marker': 'none', 'ms': 1., 'ls': '-', 'lw': 0.7, 'color': 'tab:cyan', 'label': 'IRTF/SpeX'},
                            'KPNO 4m R-C Spec': {'marker': 'none', 'ms': 1., 'ls': '-', 'lw': 0.7, 'color': 'tab:pink', 'label': 'KPNO/R-C'}}],
              xlim=(0.6, 2.6),
              ylim=(-1.2e-14, 1.2e-13),
              ylim_res=(-5., 5.),
              scale=('linear', 'linear'),
              legend=[{'loc': 'upper right', 'frameon': False, 'fontsize': 11.},
                      {'loc': (0.2, 0.1), 'frameon': False, 'fontsize': 11.}],
              figsize=(6., 2.5),
              object_type='star',
              leg_param=['teff', 'logg', 'radius'],
              units=('um', 'W m-2 um-1'),
              output='spectrum.png')

[tomasstolker]

[kelle]

This is SUPER COOL! I'm so glad you got it to work!

When you say you added the resolution manually, are you referring to spec_res = 100? Let's have this convo over in #274.

Is this new function described in the species docs somewhere?

I noticed that you wrote a helper function (_fetch_bibcode) and a bunch of other code which does "obviously" useful things. We have a companion package, astrodb_utils which can/should do a bunch of this stuff. From your perspective, is it preferable for you to use your own code or would you prefer to import astrodb_utils and use provided functions to get the bibcodes, source names, photometry, parallax, etc?

[tomasstolker]

The docstring of the function can be found here: https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database.add_simple_object

I didn't create a tutorial, but happy to add an example later on. Although for fitting spectra/photometry, this tutorial could be followed:

https://species.readthedocs.io/en/latest/tutorials/fitting_model_spectra.html

The only difference is that add_simple_object instead of add_companion/add_object is used. add_simple_object and add_object can be combined though. For example, first adding all the SIMPLE data and then adding some proprietary spectrum with add_object for the same target.

When you say you added the resolution manually, are you referring to spec_res = 100?

No that value I choose a bit arbitrarily. It is used for the model spectra that are plotted, although mainly the best-fit is seen. The random posterior spectra in gray are hardly visible since this a a high S/N spectrum and there is a somewhat small number of free parameters.

I did fit for an uncertainty inflation, to account for model systematics, but in the residuals there are clearly some systematics. For the calculation of the residuals, it uses the actual resolution of the data though, so the 100 is only used for the upper panel in the plot.

The resolution is asked as manual input with input() when a spectrum is found for the requested SIMPLE object, so it is not shown in the code above.

We have a companion package, astrodb_utils which can/should do a bunch of this stuff. From your perspective, is it preferable for you to use your own code or would you prefer to import astrodb_utils and use provided functions to get the bibcodes, source names, photometry, parallax, etc?

Using functions from astrodb_utils would certainly be fine! I simply adapted code from SEDkitSIMPLE, including the function for fetching bibcodes. If there is a function from astrodb that could be called directly then that may be more robust for any changes in SIMPLE in the future.

[kelle]

This seems to be working well! Let's open issues as we identify other ways to help these two tools work well together.