Simulating Light-Curves with a Cadence

This notebook demonstrates the simulation of normal Type Ia Supernova (SN Ia) light-curves using realistic cadences and atmospheric variability expected from LSST. To achieve this we use data from the PLAsTICC simulations to establish the cadence, light-curve parameters, and location of SNe observed by LSST. Light-curves with customized time-variable PWV transmission effects are then simulated using these parameters.

Plots in this notebook are rendered using a subset of the PLAsTICC simulation data. To access this data, configure your working environment as outlined in the docs.

import numpy as np
import sncosmo
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.table import vstack
from matplotlib import pyplot as plt
from pwv_kpno.gps_pwv import GPSReceiver

from snat_sim import models, constants as const
from snat_sim.plasticc import PLAsTICC

WARNING: AstropyDeprecationWarning: The update_default_config function is deprecated and may be removed in a future version. [sncosmo]

from snat_sim.data_paths import paths_at_init
paths_at_init.get_plasticc_dir()

PosixPath('/home/djperrefort/miniconda3/envs/SN-PWV/lib/python3.9/site-packages/data/plasticc')

plt.rcParams['figure.dpi'] = 100

The PLAsTICC Data

Instead of evaluating different cadences from scratch, we take the PLAsTICC simulations as a starting point. First we check what cadence simulations are available on the notebook’s host server.

PLAsTICC.get_available_cadences()

['alt_sched_rolling', 'alt_sched']

We load an example light-curve from one of these data sets and demonstrate the data below. Each cadence includes simulations run with multiple supernova models. In this notebook we only need simulations for normal SNe (Model 11).

alt_sched = PLAsTICC('alt_sched', 11)
demo_header_files = alt_sched.get_model_headers()
print('Available Light-curves: ', alt_sched.count_light_curves())

Available Light-curves:  115

# Data is loaded into memory in chunks. This cell may take a second.
snid, params, cadence, light_curve = next(alt_sched.iter_cadence(include_lc=True, verbose=False))

light_curve.to_astropy()

Table length=191

time	band	flux	fluxerr	zp	zpsys	phot_flag
float64	str15	float32	float32	float32	str2	int32
61954.3306	lsst_hardware_z	15.283731	20.61587	30.89	AB	0
...	...	...	...	...	...	...
62468.0277	lsst_hardware_y	6.5977235	42.13687	30.05	AB	0
62468.0439	lsst_hardware_g	1.3531196	3.310382	31.4	AB	0

sncosmo.plot_lc(light_curve.to_astropy());

Simulating Light-Curves

Since we need to add in our own atmospheric variability, the pre-tabulated flux values above are of limited use. Instead, we use the PLAsTICC meta-data to establish the cadence and model parameters for each simulated SN. This information is then used to simulate our own light-curves with sncosmo.

model_for_sim = models.SNModel('salt2-extended')
model_for_sim.update({p: v for p, v in params.items() if p in model_for_sim.param_names})
model_for_sim.set_source_peakabsmag(const.betoule_abs_mb, 'standard::b', 'AB', cosmo=const.betoule_cosmo)

cadence.zp = np.full_like(cadence.zp, 30)
duplicated_lc = model_for_sim.simulate_lc(cadence)

duplicated_lc.to_astropy()

Table length=191

time	band	flux	fluxerr	zp	zpsys	phot_flag
float64	str15	float64	float64	float64	str2	float64
61954.3306	lsst_hardware_z	54.811352066851164	47.31999969482422	30.0	AB	0.0
...	...	...	...	...	...	...
62468.0277	lsst_hardware_y	-69.66747927697624	56.18000030517578	30.0	AB	0.0
62468.0439	lsst_hardware_g	2.6161596433078604	13.300000190734863	30.0	AB	0.0

When comparing our own “duplicated” light-curves with the original PLAsTICC simulations, we care about three things: 1. Do the light-curves have the same parameters (except for x0, which depends on the spectral template and cosmology)? 1. Do the light-curves have the same time sampling? 1. Do the light-curves have a similar SNR?

The first point above is ensured by the snat_sim test suite. The second is easily verified by looking at the plotted data:

sncosmo.plot_lc(duplicated_lc.to_astropy());

When examining the SNR, it is more informative to look at a distribution over many light-curves:

def compile_light_curves(dao: PLAsTICC, override_zp: float = 30):
    """Create stacked tables for PLAsTICC simulated light-curves and their duplicated simulations

    Args:
        dao: The PLAsTICC data access object
        override_zp: Overwrite the zero-point used by plasticc with this number

    Returns:
        - Stacked table of PLAsTICC light-curves
        - Stacked table of duplicated light-curves
    """

    model_for_sim = models.SNModel('salt2-extended')
    data_iterator = dao.iter_cadence(iter_lim=1000, include_lc=True, verbose=False)

    plasticc_data = []
    duplicated_data = []
    for _, params, cadence, plasticc_lc in data_iterator:
        model_for_sim.update({p: v for p, v in params.items() if p in model_for_sim.param_names})
        model_for_sim.set_source_peakabsmag(const.betoule_abs_mb, 'standard::b', 'AB', cosmo=const.betoule_cosmo)

        if override_zp:
            cadence.zp = np.full_like(cadence.zp, override_zp)

        try:
            new_lc = model_for_sim.simulate_lc(cadence)

        except:
            continue

        else:
            plasticc_data.append(plasticc_lc.to_astropy())
            duplicated_data.append(new_lc.to_astropy())

    return vstack(plasticc_data), vstack(duplicated_data)


def plot_snr_distributions(dao, bins=np.arange(-10, 16)):
    """Plot the SNR distribution in each band for PLAsTICC light-curves and their duplicates

    Args:
        lc_list (List[Table]): List of light-curves in the PLAsTICC data model
        bins          (array): Bins to use when generating the histogram
    """

    snr_results = compile_light_curves(dao)

    fig, axes = plt.subplots(2, 5, figsize=(15, 6), sharex='col', sharey=True)
    for combined_data, axis_row in zip(snr_results, axes):
        for band, axis in zip('ugrizy', axis_row):
            band_data = combined_data[combined_data['band'] == 'lsst_hardware_' + band]
            snr = (band_data['flux'] / band_data['fluxerr'])
            avg = np.average(snr)
            std = np.std(snr)

            axis.hist(snr, alpha=.5, bins=bins)
            axis.axvline(avg, linestyle='--', color='k', label=f'Average = {avg: .2f}')
            axis.axvline(avg - std, linestyle=':', color='k', label=f'Stdev = {std: .2f}')
            axis.axvline(avg + std, linestyle=':', color='k')

            axis.set_title(band)
            axis.set_xlim(min(bins), max(bins))
            axis.legend()

    for axis in axis_row:
        axis.set_xlabel('SNR')

    axes[0, 0].set_ylabel('PLAsTICC')
    axes[1, 0].set_ylabel('SNCosmo')

    plt.show()

plot_snr_distributions(alt_sched)

Adding Atmospheric Effects

To simulate light-curves with PWV effects, we need a model for the PWV over time. This is easily defined as follows:

ctio = GPSReceiver('CTIO', data_cuts={'PWV': [(0, 25)]})
ctio.download_available_data(year=range(2012, 2018))
pwv_model = models.PWVModel.from_suominet_receiver(ctio, 2016, [2017])

Next we need a SN Model that incorporates the PWV concentration along line of sight (i.e. the PWV scaled for the airmass at the time of observation). The sncosmo package doesn’t have a clearly defined approach to adding time variable propagation effects. We use custom classes from the snat_sim package instead, which mimic sncosmo but include time variable propagation effects.

pwv_effect = models.VariablePWVTrans(pwv_model)

demo_model_with_pwv = models.SNModel(
    source='salt2-extended',
    effects=[pwv_effect],
    effect_names=[''],
    effect_frames=['obs']
)

def plot_variable_pwv_sn_model(model_with_pwv, phase=0, params=None):
    """Over plot a sncosmo model with and without temporally variable PWV

    Args:
        model_with_pwv (SNModel): sncosmo source to plot
        phase            (float): Phase of the supernova to plot
        params            (dict): Non-PWV related parameters for the model
    """

    params = params or dict()
    wave = np.arange(3000, 12000)
    time = phase + params.get('t0', 0)

    model_without_pwv = sncosmo.Model(model_with_pwv.source)
    model_without_pwv.update({k: v for k, v in params.items() if k in model_without_pwv.param_names})
    flux_without_pwv = model_without_pwv.flux(time, wave)

    model_with_pwv.update(params)
    flux_with_pwv = model_with_pwv.flux(time, wave)

    fig, ax = plt.subplots(1, 1, figsize=(6, 3))
    ax.plot(wave, flux_without_pwv, label='Base Model', color='C1')
    ax.plot(wave, flux_with_pwv, label='Model with PWV', color='C0')
    ax.set_title('Simulated Flux')
    ax.set_ylabel('Flux')
    ax.legend()
    ax.set_xlabel('Wavelength (A)')
    ax.set_xlim(min(wave), max(wave))

    plt.tight_layout()

plot_variable_pwv_sn_model(
    demo_model_with_pwv,
    params = {
        'z': 0.752069652,
        'x0': 3e-06,
        'x1': -1.8,
        'c': -0.1,
    }
)

The PWV along line of sight is equivalent to the PWV at zenith times the airmass of the observation. We quickly validate the airmass calculation of the VariablePWVTrans class.

def plot_airmass_validation(cadence, model=11, mjd=0):
    """Plot the airmass of simulated PLAsTICC SNe locations as observed from LSST

    Args:
        cadence_name (PLAsTICC): Simulated PLAsTICC cadence to plot SNe for
        mjd             (float): Date to plot airmasses for
    """

    ra = []
    dec = []

    header_path = cadence.get_model_headers()[0]
    header_data = fits.open(header_path)[1].data
    ra.extend(header_data['RA'])
    dec.extend(header_data['DECL'])

    if mjd == 'peak':
        mjd = header_data['PEAKMJD']

    airmass = models.PWVModel.calc_airmass(time=mjd, ra=ra, dec=dec, raise_below_horizon=False)
    is_positive_airmass = np.array(airmass) >= 0
    positive_airmass = airmass[is_positive_airmass]

    sn_coord = SkyCoord(ra, dec, unit=u.deg).galactic
    positive_coords = sn_coord[is_positive_airmass]
    negative_coords = sn_coord[~is_positive_airmass]

    plt.figure(figsize=(10, 5))
    plt.subplot(111, projection='aitoff')
    plt.grid(True)

    scat = plt.scatter(positive_coords.l.wrap_at('180d').radian, positive_coords.b.radian, c=positive_airmass, vmin=1, vmax=8, s=10)
    plt.scatter(negative_coords.l.wrap_at('180d').radian, negative_coords.b.radian, c='lightgrey', label='Over Horizon')
    plt.legend(framealpha=1)
    plt.colorbar(scat).set_label('Airmass', rotation=270, labelpad=15)

plot_airmass_validation(alt_sched)