Fit supernova + galaxy model on a series of spectral data cubes from the Nearby Supernova Factory.
Above: Each column shows a separate observation date for a supernova, with the last four observations occuring after the supernova faded. Top row: The data integrated over a 1000 Angstrom-wide bandpass. Bottom row: The cubefit model for the scene integrated over the same bandpass.
To install a release version:
pip install http://github.com/snfactory/cubefit/archive/v0.4.2.tar.gz
Release versions are listed here. CubeFit has the following dependencies:
- Python 2.7 or 3.3+
- scipy (for optimization)
- fitsio (for FITS file I/O)
- pyfftw (for fast FFTs)
Cubefit consists of three command-line executables.
Fit the model, write an output model to given file:
cubefit config.json output.fits
Producing subtracted data cubes is a separate step:
cubefit-subtract config.json output.fits
(This reads the input and output filenames listed in the configuration
file, and the results of the fit saved in
Finally, there are some standard plots that can be produced with
cubefit-plot config.json output.fits plotdir
plotdir is the target directory for the resulting plots).
With any command, you can run with
--help options to see all the
optional arguments, e.g.,
CubeFit requires the following keys in the input JSON file. Keys are case-sensitive. Additional keys are simply ignored.
||list||FITS data cube files|
||list||x position of MLA center relative to SN [spaxels]|
||list||y position of MLA center relative to SN [spaxels]|
||list of lists||PSF parameters for each epoch|
||list||index of final refs in lists [0-based indexing]|
||index of "master" final ref [0-based indexing]|
||list||Output subtracted cube filenames for
||list||Output SN spectrum filenames for
In each input FITS file, the following header keywords are used:
||Airmass of observation|
||Atmospheric pressure [mmHg]|
||Atmospheric temperature [degrees Celcius]|
||Parallactic angle [degrees]|
CubeFit outputs a FITS file with two extensions.
HDU 1: (image) Full galaxy model. E.g., a 32 x 32 x 779 array.
HDU 2: (binary table) Per-epoch results. Columns:
yctrFitted y position of data in model frame
xctrfitted x position of data in model frame
sn(1-d array) fitted SN spectrum in this epoch
sky(1-d array) fitted sky spectrum in this epoch
galeval(3-d array) Galaxy model evaluated on this epoch (e.g., 15 x 15 x 779)
sneval(3-d array) SN model evaluated on this epoch (e.g., 15 x 15 x 779)
The HDU 1 header contains the wavelength solution in the "CRPIX3", "CRVAL3", "CDELT3" keywords. These are simply propagated from input data cubes.
The HDU 2 header contains the fitted SN position in the model frame in the "SNX" and "SNY" kewords.
Here is an example demonstrating how to reproduce the full model of the scene for the first epoch:
>>> import fitsio >>> f = fitsio.FITS("PTF09fox_B.fits") >>> f file: PTF09fox_B.fits extension: 1 type: BINARY_TBL extname: epochs rows: 17 column info: yctr f8 xctr f8 sn f4 array sky f4 array galeval f4 array[15,15,779] sneval f4 array[15,15,779] >>> epochs = f.read() >>> i = 0 >>> scene = epochs["sky"][i, :, None, None] + epochs["galeval"][i] + epochs["sneval"][i]
What it does
A description of the original implemetation of the algorithm can be found in Bongard et al (2011).
The data are N 3-d cubes, where N are observations at different times (typically N ~ 10-30). Each cube has two spatial directions and one wavelength direction. For SNFactory, the dimensions are 15 x 15 x M, where M=779 for the blue channel M=1572 for the red channel. Each cube has both a data array and a corresponding weight (error) array.
Some of the data cubes are designated as "final refs", meaning that there is assumed to be no SN light in these epochs.
The model parameters consist of:
- 3-d galaxy model ("unconvolved"). Default size is 32 x 32 x M, where the pixel size of the model matches that of the data (both spatially and in wavelength). The model extends past the data spatially, and the spectral grid matches the data.
- N 1-d sky spectra
- N 1-d supernova spectra
- position of each of the N data cubes: 2N parameters
- position of the SN: 2 parameters
Some of these parameters are fixed. For example, the SN spectrum in each of the final refs is fixed to zero. The position of one of the final refs (the "master" final ref) is fixed to the input values because only the relative offset between data cube positions is discernable.
Additionally there are inputs that describe the point spread function (PSF) in each epoch. These are necessary for propagating the model to the data. They are derived externally and are not varied in the fit. The PSF is "3-d": it is defined on a grid of wavelengths. The spatial shape can vary between wavelengths, as can the position of the center of the PSF. The change in the center position with wavelength encodes the amount of atmospheric differential refraction (ADR).
Finally, there are two regularization parameters that determine the penalty on the galaxy model for being "rougher" (having a high pixel-to-pixel variation).
These steps are carried out in
cubefit.main.cubefit() (which is called from
the command-line script).
- Read configuration file
- Read in datacubes
- Initialize the PSF model for each epoch, based on values in the configuration file and image header.
- Initialize all model parameters to zero, except sky. For the sky, we make an initial heuristic guess based on only the data.
- Initialize the regularization, based on the input regularization parameters and a rough guess at the average galaxy spectrum.
Fit the galaxy model to just the master final ref
The position of the master final ref is fixed with respect to the model, so there is no position fit here. Also, the sky of the master final ref is fixed to the initial guess (as the galaxy model and sky level are degenerate).
Fit the position and sky of the remaining final refs
Re-fit galaxy model using all final refs
Fit position, sky and SN position of all the non-final refs
Internal API documentation
This is intended as a rough guide to the major components of the code.
||Entry point for
||Entry point for
||Entry point for
|Data structure and I/O|
||Container for data and weight arrays|
||Read a two-HDU FITS file into a DataCube|
||A 3-d PSF model made up of a Gaussian + Mofffat profile.|
||A 3-d PSF model defined by a 3-d array.|
||Instatiate a 3-d PSF model based on SNFactory-specific parameterization.|
||Evaluate a 3-d Gaussian + Moffat profile on a 3-d array.|
ADR class from the SNfactory Toolbox package and the
Hyper_PSF3D_PL class from libExtractStar are used in
||Callable that returns the penalty and gradient on it.|
||Guess sky based on lower signal spaxels compatible with variance|
||Fit the galaxy model to a single epoch of data.|
||Fit the galaxy model to multiple data cubes.|
||Fit data position and sky for a single epoch (fixed galaxy model).|
||Fit data pointing (nepochs), SN position (in model frame), SN amplitude (nepochs), and sky level (nepochs).|
||Return a figure showing data and model for all epochs.|
||Make a more detailed figure for a single epoch.|
These functions are called by
If you've clone the repository (rather than using pip to install), you
can run tests with
setup.py test. Requires the
(available via pip or conda).
All code in this repository is released under the MIT license. However, since CubeFit uses FFTW, which is GPL-licensed, the software as a whole is bound by the terms of the GPLv2.
Practically, this means that one can copy the code here, remove the dependence on FFTW and release the code under the MIT license (or other permissive license).