Generative Fitting Pipeline¶
Introduction¶
This module includes tools to generate synthetic white dwarf spectra using theoretical models and a speedy generative neural network.
These can be used to fit observed spectra in an MCMC framework to recover parameters like effective temperature, surface gravity, and radial velocity.
The synthetic spectra are generated from models developed by [2], and we use the emcee
fitting routine from [1].
Say we have an un-normalized DA spectrum from the Sloan Digital Sky Survey (SDSS), defined by wavelengths wl
, fluxes fl
, and an inverse variance array ivar
.
import wdtools
gfp = wdtools.GFP(resolution = 3) # Spectral dispersion in Angstroms
labels, e_labels, redchi = gfp.fit_spectrum(wl, flux, ivar,
mcmc = True, nwalkers = 50, burn = 100, ndraws = 100,
make_plot = True, plot_corner = True
)
The GFP first estimates the radial velocity of the provided spectrum using the H-alpha absorption line. The synthetic models are shifted to this RV during the fitting process–hence there is no interpolation or re-binning of the observed spectrum, preventing correlated errors. The spectrum is then spline-normalized, during which the strong Balmer lines are masked out. By default, only the Balmer lines are used in the likelihood. You can select which Balmer lines to include in the fit with the lines
argument (defaults to all lines from ‘alpha’ to ‘h8’).
Then, the lmfit
least squares routine is used to fit the stellar parameters. The fit is initialized at several equidistant initial temperatures governed by the nteff
keyword. The fit with the lowest chi-square is selected and returned.
If the polyorder
argument is greater than zero, then the continuum-normalized synthetic spectra have a Chebyshev polynomial of order polyorder
added to them during the fitting process. These coefficients are also solved for and returned by the GFP. If you use this option, it’s recommended to also run mcmc
so that these coefficients are properly marginalized over. Also, it’s only recommended to use polyorder
if fullspec
is True
.
If the mcmc
argument is True
, then the lmfit
estimates are used as a starting point for a full MCMC run, governed by the nwalkers
, burn
, and ndraws
arguements. If mcmc
is False
, then the returned uncertainties are estimated from the covariance matrix returned by lmfit
. Turning off mcmc
results in much quicker results, but the error estimates might not be robust.
label
and e_labels
respectively are arrays of the fitted parameters. The radial velocity (and RV error) are always the last elements in the array, so if polyorder
> 0, the label array will have temperature, surface gravity, the Chebyshev coefficients, and then RV. redchi
is the reduced chi-square statistic, which can be used as a rough estimate of the goodness-of-fit. More details are in our paper, and a more complete example can be found in the tutorial.
- 1
Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman. emcee: The MCMC Hammer. PASP, 125(925):306, March 2013. doi:10.1086/670067.
- 2
D. Koester. White dwarf spectra and atmosphere models. Mem. SAI, 81:921–931, 2010.
API¶
- class wdtools.GFP(resolution=3, specclass='DA')
Generative Fitting Pipeline.
- fit_spectrum(wl, fl, ivar=None, prior_teff=None, mcmc=False, fullspec=False, polyorder=0, norm_kw={'k': 1, 'niter': 0, 'sfac': 0.5}, nwalkers=25, burn=25, ndraws=25, threads=1, progress=True, plot_init=False, make_plot=True, plot_corner=False, plot_corner_full=False, plot_trace=False, savename=None, DA=True, crop=(3600, 7500), verbose=True, lines=['alpha', 'beta', 'gamma', 'delta', 'eps', 'h8'], lmfit_kw={'epsfcn': 0.1, 'method': 'leastsq'}, rv_kw={'distance': 100, 'edge': 15, 'nmodel': 2, 'plot': False}, nteff=3, rv_line='alpha', corr_3d=False)
Main fitting routine, takes a continuum-normalized spectrum and fits it with MCMC to recover steller labels.
- Parameters
wl (array) – Array of observed spectral wavelengths
fl (array) – Array of observed spectral fluxes, continuum-normalized. We recommend using the included normalize_balmer function from wdtools.spectrum to normalize DA spectra, and the generic continuum_normalize function for DB spectra.
ivar (array) – Array of observed inverse-variance for uncertainty estimation. If this is not available, use ivar = None to infer a constant inverse variance mask using a second-order beta-sigma algorithm. In this case, since the errors are approximated, the chi-square likelihood may be inexact - treat returned uncertainties with caution.
prior_teff (tuple, optional) – Tuple of (mean, sigma) to define a Gaussian prior on the effective temperature parameter. This is especially useful if there is strong prior knowledge of temperature from photometry. If not provided, a flat prior is used.
mcmc (bool, optional) – Whether to run MCMC, or simply return the errors estimated by LMFIT
fullspec (bool, optional) – Whether to fit the entire continuum-normalized spectrum, or only the Balmer lines.
polyorder (int, optional) – Order of additive Chebyshev polynomial during the fitting process. Can usually leave this to zero unless the normalization is really bad.
norm_kw (dict, optional) – Dictionary of keyword arguments that are passed to the spline normalization routine.
nwalkers (int, optional) – Number of independent MCMC ‘walkers’ that will explore the parameter space
burn (int, optional) – Number of steps to run and discard at the start of sampling to ‘burn-in’ the posterior parameter distribution. If intitializing from a high-probability point, keep this value high to avoid under-estimating uncertainties.
ndraws (int, optional) – Number of ‘production’ steps after the burn-in. The final number of posterior samples will be nwalkers * ndraws.
threads (int, optional) – Number of threads for distributed sampling.
progress (bool, optional) – Whether to show a progress bar during the MCMC sampling.
plot_init (bool, optional) – Whether to plot the continuum-normalization routine
make_plot (bool, optional) – If True, produces a plot of the best-fit synthetic spectrum over the observed spectrum.
plot_corner (bool, optional) – Makes a corner plot of the fitted stellar labels
plot_corner_full (bool, optional) – Makes a corner plot of all sampled parameters, the stellar labels plus any Chebyshev coefficients if polyorder > 0
plot_trace (bool, optiomal) – If True, plots the trace of posterior samples of each parameter for the production steps. Can be used to visually determine the quality of mixing of the chains, and ascertain if a longer burn-in is required.
savename (str, optional) – If provided, the corner plot and best-fit plot will be saved as PDFs in the working folder.
DA (bool, optional) – Whether the star is a DA white dwarf or not. As of now, this must be set to True.
crop (tuple, optional) – The region to crop the supplied spectrum before proceeding with the fit. Can be used to exclude low-SN regions at the edge of the spectrum.
verbose (bool, optional) – If True, the routine prints several progress statements to the terminal.
lines (array, optional) – List of Balmer lines to utilize in the fit. Defaults to all from H-alpha to H8.
lmfit_kw (dict, optional) – Dictionary of keyword arguments to the LMFIT solver
rv_kw (dict, optional) – Dictionary of keyword arguments to the RV fitting routine
nteff (int, optional) – Number of equidistant temperatures to try as initialization points for the minimization routine.
rv_line (str, optional) – Which Balmer line to use for the radial velocity fit. We recommend ‘alpha’.
corr_3d (bool, optional) – If True, applies 3D corrections from Tremblay et al. (2013) to stellar parameters before returning them.
- Returns
Returns the fitted stellar labels along with a reduced chi-square statistic with the format: [[labels], [e_labels], redchi]. If polyorder > 0, then the returned arrays include the Chebyshev coefficients. The radial velocity (and RV error) are always the last elements in the array, so if polyorder > 0, the label array will have temperature, surface gravity, the Chebyshev coefficients, and then RV.
- Return type
array
- inv_label_sc(label_array)
Inverse label scaler to transform Teff and logg from [0,1] to original scale based on preset bounds.
- Parameters
label_array (array) – Scaled array with Teff in the first column and logg in the second column
- Returns
Unscaled array
- Return type
array
- label_sc(label_array)
Label scaler to transform Teff and logg to [0,1] interval based on preset bounds.
- Parameters
label_array (array) – Unscaled array with Teff in the first column and logg in the second column
- Returns
Scaled array
- Return type
array
- spectrum_sampler(wl, teff, logg, *polyargs, specclass=None)
Wrapper function that talks to the generative neural network in scaled units, and also performs the Gaussian convolution to instrument resolution.
- Parameters
wl (array) – Array of spectral wavelengths on which to generate the synthetic spectrum
teff (float) – Effective surface temperature of sampled spectrum
logg (float) – log surface gravity of sampled spectrum (cgs)
polyargs (float, optional) – All subsequent positional arguments are assumed to be coefficients for the additive Chebyshev polynomial. If none are provided, no polynomial is added to the model spectrum.
specclass (str, optional) – Whether to use hydrogen-rich (DA) or helium-rich (DB) atmospheric models. If none, reverts to default.
- Returns
Synthetic spectrum with desired parameters, interpolated onto the supplied wavelength grid and convolved with the instrument resolution.
- Return type
array
- spline_norm_DA(wl, fl, ivar, kwargs={'k': 3, 'niter': 3, 'sfac': 1}, crop=None)
Masks out Balmer lines, fits a smoothing spline to the continuum, and returns a continuum-normalized spectrum
- Parameters
wl (array) – Array of observed spectral wavelengths.
fl (array) – Array of observed spectral fluxes.
ivar (array) – Array of observed inverse-variance.
kwargs (dict, optional) – Keyword arguments that are passed to the spline normalization function
crop (tuple, optional) – Defines a start and end wavelength to crop the spectrum to before continuum-normalization.
- Returns
If crop is None, returns a 2-tuple of (normalized_flux, normalized_ivar). If a crop region is provided, then returns a 3-tuple of (cropped_wavelength, cropped_normalized_flux, cropped_normalized_ivar).
- Return type
tuple
- synth_spectrum_sampler(wl, teff, logg, rv, specclass=None)
Generates synthetic spectra from labels using the neural network, translated by some radial velocity. These are _not_ interpolated onto the requested wavelength grid; The interpolation is performed only one time after the Gaussian convolution with the instrument resolution in GFP.spectrum_sampler. Use GFP.spectrum_sampler in most cases.
- Parameters
wl (array) – Array of spectral wavelengths (included for completeness, not used by this function)
teff (float) – Effective surface temperature of sampled spectrum
logg (float) – log surface gravity of sampled spectrum (cgs)
rv (float) – Radial velocity (redshift) of sampled spectrum in km/s
specclass (str ['DA', 'DB']) – Whether to use hydrogen-rich (DA) or helium-rich (DB) atmospheric models. If None, uses default.
- Returns
Synthetic spectrum with desired parameters, interpolated onto the supplied wavelength grid.
- Return type
array