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