Using a Gaussian Process model object

Introduction

This notebook will demonstrate several features of the SNModel class, which stores the model generated from the multidimensional Gaussian Process regression over a collection of transients.

These models by default are saved as .fits files and contain several important data products that will be useful to users interested in working with them:

1. SNModel.surface – the model, or “template”, of the transient collection. A prediction of the time-evolution of the class’s spectral energy distribution at an arbitrary wavelength and phase. It is either an sklearn GaussianProcessRegressor object or a numpy array containing the median SED surface and its IQR.

2. template_grid – a grid of photometry used to normalize the photometry of the collection of objects fit by the GP. This is used to calculate the residuals of the photometry as input to the GPs, and is needed when converting between “observed” photometry and photometric residuals used by the GP.

3. phase_grid and wavelength_grid – grids in wavelength and phase space, covering the range of values fit by the GP, with shapes corresponding to the output shape of the surface and template_grid objects.

4. sn or sn_collection – the collection of objects fit by the GP

5. norm_set – the collection of objects used to normalize the photometry

Setup

Let’s import the necessary objects.

from caat import SNModel, SN
import numpy as np

import logging

logger = logging.getLogger()
logger.setLevel(logging.ERROR)

Configuration file not found:

    /Users/craigpellegrino/.dustmapsrc

To create a new configuration file in the default location, run the following python code:

    from dustmaps.config import config
    config.reset()

Note that this will delete your configuration! For example, if you have specified a data directory, then dustmaps will forget about its location.

1. Load a saved model object

First let’s load a saved GP template object. If a file already exists for an object, providing the filename will automatically load in the saved model object and all other necessary data products.

model = SNModel(
    surface="SESNe_SNIIb_GP_model.fits",
)

The SNModel object has a number of useful methods, which we’ll demonstrate now.

2. Predict a light curve and compare a light curve to observed photometry

One straightforward usage of an SNModel object is to use the Gaussian Process regression to predict a light curve. We’ll give the object a min and max phase as well as a wavelength (in Angstroms):

model.predict_lightcurve(-20.0, 45.0, 5000)

The above plot shows a predicted light curve (dark blue line) and the uncertainty in the prediction (light blue lines). We can also see how this prediction stacks up against actual photometry for an object in this class:

model.predict_lightcurve(-20.0, 45.0, 5000, show=False)
gkg = SN(name="SN2016gkg")
model.compare_lightcurve_with_photometry(gkg, filt="V")

Not bad! This is one of the more peculiar objects in this sample, SN 2016gkg, which has a strong early-time light curve bump due to shock cooling emission. Even so, the photometric evolution is captured within the uncertainty of the GP fairly well, and the rise to the main light curve peak is reproduced.

3. Predict an SED and compare an SED to observed spectra

Another basic usage is to predict an SED using a GP model at an arbitrary phase, as this next example shows:

model.predict_sed(
    wavelength_min=2500,
    wavelength_max=8000,
    phase=0.0,
    show=True
)

If desired, you can also compare the SED predictions to low resolution spectra. However, this isn’t an exact comparison since the spectra may not be flux calibrated, so the conversion from flux units to normalized flux isn’t exact, but it gives us a rough estimate if the shape of the SED matches observations.

4. Produce synthetic photometry for a class of objects

These GP models can also be used to produce synthetic photometry of a given transient class at an arbitrary phase and wavelength. This is especially useful for constructing training sets at high redshifts to better classify transients discovered by Rubin LSST and Roman.

model = SNModel(
    surface="SESNe_SNIIb_GP_model.fits",
)

model.predict_photometry_points(
    wavelengths = np.asarray([5000.0, 5000.0, 5000.0, 7000.0, 7000.0, 7000.0]),
    phases = np.asarray([-5.0, 0.0, 10.0, -5.0, 0.0, 10.0]),
    show=True
)

(array([-5.,  0., 10., -5.,  0., 10.]),
 array([-0.06825227,  0.38219782, -1.8790352 ,  0.46723262,  1.00798788,
         0.05463794]),
 array([0.41772546, 0.32090284, 0.42934931, 0.58340145, 0.47485424,
        0.58166458]))

5. Fit input photometry of a new object

Maybe you want to see how well this model performs for new photometry. You can input a dictionary or pandas dataframe of photometry and plot a fit, even extrapolated out to phases beyond those covered by the input. The input photometry magnitudes should be relative to the light curve peak–in the example below, the first value with magnitude = -1.0 means that point is 1 magnitude fainter than the peak magnitude.

mock_photometry_values = {
    "Filter": ["B", "B", "B", "B", "V", "V", "V", "V", "i", "i", "i", "i"],
    "Phase": [-15.0, -10.0, -4.0, 1.0, -15.0, -10.0, -4.0, 1.0, -15.0, -10.0, -4.0, 1.0],
    "Mag": [-1.0, -0.6, -0.3, 0.3, -1.1, -0.7, -0.25, 0.1, -1.2, -0.5, -0.35, -0.1],
    "MagErr": [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
}

model.fit_photometry(
    photometry=mock_photometry_values,
    show=True,
    phase_max=1.0,
)

Let’s run this again, but this time extrapolate the fit out and see how it performs.

model.fit_photometry(
    photometry=mock_photometry_values,
    phase_min=-15,
    phase_max=45,
    show=True
)

Or maybe it would be useful to visualize the uncertainty as a series of different light curves, rather than an uncertainty region. We can plot multiple individual fits per filter using the nsamples parameter:

model.fit_photometry(
    photometry=mock_photometry_values,
    phase_min=-15,
    phase_max=45,
    show=True,
    nsamples=10
)

6. Fit photometry of a SN object

The SNModel objects can also be used to produce fits to existing SN objects. The SNModel.fit_photometry method can also take a SN object as an argument via the sn_to_fit parameter. In this case, it’ll default to fitting the passed-in SN object, rather than any optionally-provided photometry. Users have the ability to specify what filters to fit, or to use the default setting of fitting each filter for which there is photometry. For example, let’s use the SN IIb template model to fit the g- and r-band photometry of a double-peaked SN IIb:

sn = SN(name="SN2020rsc")

model.fit_photometry(
    sn_to_fit=sn,
    filters_to_fit=["g", "r"],
    show=True,
    nsamples=1
)

7. Predict photometry for a SN outside its data range

Finally, SNModel can be used to generate new synthetic photometry that is outside a SN object’s coverage in phase or wavelength. This is done by taking advantage of the fit_photometry method. This method has a flag keep_new_fit which, if set to True, uses the generated GPR fit to the input photometry for all future method calls of this class instance.

Let’s see an example of this. We’ll load the SNModel object for Type IIb SNe, then fit the Type IIb SN 2020rsc again. This object only has data in the optical (g- through i-band) is missing data in the Swift filters; however, using the SNModel object fit to the object’s photometry, we can predict a missing filter’s photometry:

from caat.utils import WLE


sn = SN(name="SN2020rsc")

model.fit_photometry(
    sn_to_fit=sn,
    show=True,
    nsamples=1,
    keep_new_fit=True
)

model.predict_lightcurve(
    phase_min=-18,
    phase_max=45,
    wavelength=WLE["U"],
    show=True,
)

The new fit has extrapolated a light curve in a new band! As you can see, the shape closely follows the shapes of the light curves in the other filters, with reasonably large error bars at late phases. Similarly, you can use the predict_photometry_points method in a similar manner. This demonstrates the predictive power of these GPR models!

Next steps:

These GP model surfaces can be manipulated in a number of ways. The SNModel object saves all information needed to generate the underlying sklearn GaussianProcessRegressor object, which allows for repeated GP fits to new or existing data. The above features are just a sample of what can be done with these novel models.