Basic SN Ia Example

In this notebook we look at a demo simulation and plotting of a Type Ia supernova. This notebook is intended to serve as an example of how users can set up a model and interact with the results.

import numpy as np
import pandas as pd

from lightcurvelynx.astro_utils.passbands import PassbandGroup
from lightcurvelynx.astro_utils.snia_utils import DistModFromRedshift, X0FromDistMod
from lightcurvelynx.graph_state import GraphState
from lightcurvelynx.math_nodes.np_random import NumpyRandomFunc
from lightcurvelynx.models.sncosmo_models import SncosmoWrapperModel
from lightcurvelynx.obstable.lsst_obstable import LSSTObsTable
from lightcurvelynx.simulate import compute_single_noise_free_lightcurve, simulate_lightcurves
from lightcurvelynx.survey_info import SurveyInfo
from lightcurvelynx.utils.extrapolate import LinearDecay
from lightcurvelynx.utils.plotting import plot_lightcurves

/home/docs/checkouts/readthedocs.org/user_builds/lightcurvelynx/envs/stable/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Fake Survey Data

For the purposes of this notebook, we create a very basic survey using the LSST instrument that looks at one region of the sky every night in a single filter. While this is not a realistic cadence, it provides a useful baseline for examining the behavior of a type Ia supernova.

For more realistic surveys, see the examples for Rubin OpSim, Rubin CCDVisit Table, or Custom Surveys.

filters = ["g", "r", "i", "z"]
num_days = 1000

times = np.arange(num_days) + 60676.0  # MJD starting at 60676.0
fake_data = {
    "expMidptMJD": times,  # MJD
    "ra": np.full(num_days, 150.0),  # Degrees
    "dec": np.full(num_days, -30.0),  # Degrees
    "band": [filters[i % len(filters)] for i in range(num_days)],
    "exptime": np.full(num_days, 30.0),  # Seconds
    # Noise information is random from distributions around DP2 values.
    "pixelScale": np.full(num_days, 0.2),  # Arcseconds per pixel
    "seeing": np.random.uniform(low=0.84, high=0.88, size=num_days),  # arcseconds
    "skyBg": np.random.uniform(low=1420.0, high=1430.0, size=num_days),  # ADU
    "zeroPoint": np.random.normal(loc=30.0, scale=0.5, size=num_days),  # magnitudes to produce 1 count (ADU)
}
obs_table = LSSTObsTable(pd.DataFrame(fake_data))
obs_table.head()

	time	ra	dec	filter	exptime	pixel_scale	seeing	sky_bg_adu	zp_mag_adu	zp	psf_footprint	sky_bg_e
0	60676.0	150.0	-30.0	g	30.0	0.2	0.866456	1427.875380	29.682146	3.050569	42.533143	2277.461231
1	60677.0	150.0	-30.0	r	30.0	0.2	0.842861	1420.881920	29.979774	2.319154	40.248211	2266.306663
2	60678.0	150.0	-30.0	i	30.0	0.2	0.875774	1421.079104	30.673931	1.223678	43.452910	2266.621171
3	60679.0	150.0	-30.0	z	30.0	0.2	0.874234	1424.646493	29.471214	3.704707	43.300165	2272.311156
4	60680.0	150.0	-30.0	g	30.0	0.2	0.854505	1427.627125	29.742136	2.886589	41.367911	2277.065264

Filter Information

For this example, we use the passbands from the LSST preset (but use a cached local version in the test directory to avoid a download). In most cases users will want to use data/passbands/ from the root directory.

# Use a (possibly older) cached version of the passbands to avoid downloading them.
table_dir = "../../tests/lightcurvelynx/data/passbands"
passband_group = PassbandGroup.from_preset(preset="LSST", table_dir=table_dir)
passband_group.plot()

The combination of the observations table and the bandpass filter provide all the survey information we need for the simulation. Next up is the model itself.

Create the model

To generate simulated light curves we need to define the properties of the object from which to sample. In this notebook, we use sncosmo’s SALT2 model for Type Ia supernova.

# Sample redshift uniformly between 0.01 and 0.2.
redshift_sampler = NumpyRandomFunc("uniform", low=0.01, high=0.2)

# Use given values the cosmological parameters (H0=73.0, Omega_m=0.3).
# Then compue the distance modulus from the redshift (taking the redshift sampler as input).
distmod_func = DistModFromRedshift(redshift_sampler, H0=73.0, Omega_m=0.3)

# Sample x1, c, and m_abs from distributions motivated by typical SNIa values.
x1_func = NumpyRandomFunc("normal", loc=0, scale=2.0)
c_func = NumpyRandomFunc("normal", loc=0, scale=0.02)
m_abs_func = NumpyRandomFunc("normal", loc=-19.3, scale=0.1)

# Compute x0 from the other parameters using the standard Tripp formula,
# and the distance modulus from above.
x0_func = X0FromDistMod(
    distmod=distmod_func,  # Use the computed distance modulus from redshift as input.
    x1=x1_func,  # Use the sampled x1 values as input.
    c=c_func,  # Use the sampled c values as input.
    alpha=0.14,  # Use a constant alpha value motivated by typical SNIa fits.
    beta=3.1,  # Use a constant beta value motivated by typical SNIa fits.
    m_abs=m_abs_func,  # Use the sampled m_abs values as input.
)

# t0 for the super nova is sampled uniformly over the time range of the observations.
t0_func = NumpyRandomFunc("uniform", low=np.min(times), high=np.max(times))

Sncomso’s SALT2 model is only defined over a limited range of phase space (times). So we need to decide how to handle points outside that range. We don’t want the flux to immediately drop to zero, because that produces unrealistic artifacts. Instead we use one of the extrapolation functions built into LightCurveLynx. In this case we will use a linear decay to zero over 50 days.

More types of extrapolation and more discussion of how to use them can be found in the extrapolation notebook.

time_extrapolation = LinearDecay(50.0)

We then define the model using the SncosmoWrapperModel class. All of the parameters are set using the samplers defined above. We use a constant (RA, Dec) so that the object is visible in every image. For most simulations, you will want to draw the position of the object randomly from the survey footprint. See the position sampling notebook for examples.

source = SncosmoWrapperModel(
    "salt2-h17",  # Model name
    t0=t0_func,
    x0=x0_func,
    x1=x1_func,
    c=c_func,
    ra=150.0,
    dec=-30.0,
    redshift=redshift_sampler,
    node_label="source",
    time_extrapolation=time_extrapolation,
)

Generate the simulations

We can now generate random simulations with all the information defined above. The light curves are written in the nested-pandas format for easy analysis.

survey_info = SurveyInfo(obstable=obs_table, passbands=passband_group, survey_name="LSST")
lightcurves = simulate_lightcurves(
    source,  # The model we are simulating.
    5,  # The number of simulations to run,
    survey_info,  # The survey info object containing the ObsTable and PassbandGroup.
    rest_time_window_offset=(-50, 100),  # Simulate from -50 to +100 days of t0 in the rest frame.
)

# Drop the parameters column to make the output more concise, but save it for generating
# the noise-free lightcurves in the plotting section below.
params_col = lightcurves["params"].copy()
lightcurves = lightcurves.drop(columns=["params"])
lightcurves.head()

Simulating: 100%|██████████| 5/5 [00:00<00:00, 290.65obj/s]

	id	ra	dec	nobs	t0	z	lightcurve
0	0	150.000000	-30.000000	164	61080.100048	0.097261	mjd filter ... obs_idx is_saturated 61026.0 i ... 350 False +163 rows ... ... ... ...
1	1	150.000000	-30.000000	171	61427.309676	0.146109	mjd filter ... obs_idx is_saturated 61371.0 z ... 695 False +170 rows ... ... ... ...
2	2	150.000000	-30.000000	176	61006.093213	0.172372	mjd filter ... obs_idx is_saturated 60948.0 g ... 272 False +175 rows ... ... ... ...
3	3	150.000000	-30.000000	100	61633.435621	0.151396	mjd filter ... obs_idx is_saturated 61576.0 g ... 900 False +99 rows ... ... ... ...
4	4	150.000000	-30.000000	168	60797.792950	0.118436	mjd filter ... obs_idx is_saturated 60742.0 i ... 66 False +167 rows ... ... ... ...

5 rows x 7 columns

Now let’s plot each light curve. For reference we will also plot the underlying noise-free curve.

for idx in range(5):
    # Extract the row for this object.
    lc = lightcurves.iloc[idx]

    # Unpack the nested columns (filters, mjd, flux, and flux error). This is the
    # data from the simulation itself.
    lc_filters = np.asarray(lc["lightcurve"]["filter"], dtype=str)
    lc_mjd = np.asarray(lc["lightcurve"]["mjd"], dtype=float)
    lc_flux = np.asarray(lc["lightcurve"]["flux_perfect"], dtype=float)
    lc_fluxerr = np.asarray(lc["lightcurve"]["fluxerr"], dtype=float)

    # Extract the parameters used during the simulation of this object (stored in the "params" column).
    # Use that to compute the noise-free light curve for this object, which we will plot alongside the
    # simulated light curve.
    params_dict = params_col.iloc[idx]
    noise_free_lcs = compute_single_noise_free_lightcurve(
        source,
        GraphState.from_dict(params_dict),
        passband_group,
        rest_frame_phase_min=-50.0,  # 50 days before t0
        rest_frame_phase_max=100.0,  # 100 days after t0
        rest_frame_phase_step=0.5,  # 2 samples per day
    )

    # Extract some parameters for the title of the plot.
    redshift = params_dict["source.redshift"]
    t0 = params_dict["source.t0"]

    plot_lightcurves(
        fluxes=lc_flux,
        times=lc_mjd,
        fluxerrs=lc_fluxerr,
        filters=lc_filters,
        underlying_model=noise_free_lcs,
        title=f"Simulated SNIa Light Curve (redshift={redshift:.4f}, t0={t0:.1f})",
    )