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.
[1]:
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.
[2]:
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()
[2]:
| 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.872924 | 1426.199882 | 29.634637 | 3.187018 | 43.170501 | 2274.788811 |
| 1 | 60677.0 | 150.0 | -30.0 | r | 30.0 | 0.2 | 0.861598 | 1421.463347 | 29.641733 | 3.166258 | 42.057486 | 2267.234038 |
| 2 | 60678.0 | 150.0 | -30.0 | i | 30.0 | 0.2 | 0.859972 | 1425.779463 | 30.209083 | 1.877613 | 41.898952 | 2274.118243 |
| 3 | 60679.0 | 150.0 | -30.0 | z | 30.0 | 0.2 | 0.841720 | 1425.670139 | 29.018030 | 5.623766 | 40.139328 | 2273.943871 |
| 4 | 60680.0 | 150.0 | -30.0 | g | 30.0 | 0.2 | 0.859803 | 1421.547370 | 30.381781 | 1.601500 | 41.882502 | 2267.368055 |
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.
[3]:
# 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.
[4]:
# 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.
[5]:
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.
[6]:
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,
)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/lightcurvelynx/envs/stable/lib/python3.12/site-packages/lightcurvelynx/models/sncosmo_models.py:78, in SncosmoWrapperModel.__init__(self, source_name, node_label, wave_extrapolation, time_extrapolation, seed, **kwargs)
77 try:
---> 78 from sncosmo.models import get_source
79 except ImportError as err: # pragma: no cover
ModuleNotFoundError: No module named 'sncosmo'
The above exception was the direct cause of the following exception:
ImportError Traceback (most recent call last)
Cell In[6], line 1
----> 1 source = SncosmoWrapperModel(
2 "salt2-h17", # Model name
3 t0=t0_func,
4 x0=x0_func,
File ~/checkouts/readthedocs.org/user_builds/lightcurvelynx/envs/stable/lib/python3.12/site-packages/citation_compass/citation.py:73, in CiteClass.__init_subclass__.<locals>.init_wrapper(*args, **kwargs)
69 @wraps(original_init)
70 def init_wrapper(*args, **kwargs):
71 # Save the citation as USED when it is first called.
72 CITATION_COMPASS_REGISTRY.mark_used(cls._citation_compass_name)
---> 73 return original_init(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/lightcurvelynx/envs/stable/lib/python3.12/site-packages/lightcurvelynx/models/sncosmo_models.py:80, in SncosmoWrapperModel.__init__(self, source_name, node_label, wave_extrapolation, time_extrapolation, seed, **kwargs)
78 from sncosmo.models import get_source
79 except ImportError as err: # pragma: no cover
---> 80 raise ImportError(
81 "sncosmo package is not installed by default. To use the SncosmoWrapperModel, "
82 "please install sncosmo. For example, you can install it with "
83 "`pip install sncosmo` or `conda install conda-forge::sncosmo`."
84 ) from err
86 # We explicitly ask for and pass along the PhysicalModel parameters such
87 # as node_label and wave_extrapolation so they do not go into kwargs
88 # and get added to the sncosmo model below.
89 super().__init__(
90 node_label=node_label,
91 wave_extrapolation=wave_extrapolation,
(...) 94 **kwargs,
95 )
ImportError: sncosmo package is not installed by default. To use the SncosmoWrapperModel, please install sncosmo. For example, you can install it with `pip install sncosmo` or `conda install conda-forge::sncosmo`.
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.
[7]:
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()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 3
1 survey_info = SurveyInfo(obstable=obs_table, passbands=passband_group, survey_name="LSST")
2 lightcurves = simulate_lightcurves(
----> 3 source, # The model we are simulating.
4 5, # The number of simulations to run,
5 survey_info, # The survey info object containing the ObsTable and PassbandGroup.
6 rest_time_window_offset=(-50, 100), # Simulate from -50 to +100 days of t0 in the rest frame.
NameError: name 'source' is not defined
Now let’s plot each light curve. For reference we will also plot the underlying noise-free curve.
[8]:
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})",
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 3
1 for idx in range(5):
2 # Extract the row for this object.
----> 3 lc = lightcurves.iloc[idx]
4
5 # Unpack the nested columns (filters, mjd, flux, and flux error). This is the
6 # data from the simulation itself.
NameError: name 'lightcurves' is not defined