Building Simple Models

In this tutorial we look at how to build a simple model and sample the parameters from a variety of basic sources. For more complex sampling see the notebooks: sampling.ipynb, adding_models.ipynb, and using_pzflow.ipynb.

Parameterized Nodes

All sources of information in LightCurveLynx live as ParameterizedNode objects. This allows us to link the nodes (and their variables) together and sample them as a single block. As you will see in this tutorial, most of the nodes are specific to the object that you want to simulate. For example if we wanted to create a model with constant brightness in the night sky with a brightness of 10, we could use:

import matplotlib.pyplot as plt
import numpy as np

from lightcurvelynx.models.basic_models import ConstantSEDModel

model = ConstantSEDModel(brightness=10.0, node_label="a_star")

/home/docs/checkouts/readthedocs.org/user_builds/lightcurvelynx/envs/latest/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

Users can get information about the parameterized node using Python’s builtin help() function or display a list of the parameters with the describe_params() function.

model.describe_params()

Parameters in a_star:

ra:
    Description: The object's right ascension (in degrees)
    Source: CONSTANT with value = None

dec:
    Description: The object's declination (in degrees)
    Source: CONSTANT with value = None

redshift:
    Description: The object's redshift.
    Source: CONSTANT with value = None

t0:
    Description: The phase offset in MJD.
    Source: CONSTANT with value = None

distance:
    Description: The object's luminosity distance (in pc)
    Source: CONSTANT with value = None

brightness:
    Description: The inherent brightness of the object.
    Source: CONSTANT with value = 10.0

A ParameterizedNode can then be sampled with the sample_parameters() function. This samples the internal parameters of the model. For example we might be sampling an object’s brightness (as with the constant SED model), a host galaxy’s mass, or even the (RA, dec) position of the observations. The sample_parameters() function returns a special data structure, the GraphState, that captures this state and can later be fed into the models to generate fluxes.

Note: Users do not need to know the details of the GraphState storage, only that it can be accessed like a dictionary using the node’s label and the variable name.

state = model.sample_parameters(num_samples=10)
state["a_star"]["brightness"]

array([10., 10., 10., 10., 10., 10., 10., 10., 10., 10.])

The sample function produced 10 independent samples of our system’s state.

The brightness values of these samples are not particularly interesting because we were sampling from a fixed parameter. The brightness is always 10.0. However LightCurveLynx allows the user to set a node’s parameter from a variety of sources including constants (as with 10.0 above), the values stored in other nodes, or even the results of a “functional” or “computation” type node (more about that later).

LightCurveLynx includes the built-in NumpyRandomFunc which will sample from a given numpy function and use the results to set a given parameter.

from lightcurvelynx.math_nodes.np_random import NumpyRandomFunc

brightness_func = NumpyRandomFunc("uniform", low=11.0, high=15.5)
model2 = ConstantSEDModel(brightness=brightness_func, node_label="a_star_2")
state = model2.sample_parameters(num_samples=10)

state["a_star_2"]["brightness"]

array([11.78562231, 11.79757401, 11.49950331, 14.93927065, 12.89861205,
       13.57531371, 15.21599639, 15.3566499 , 14.10435465, 13.08463892])

Now each of our 10 samples use different a different brightness value.

We can make the distributions of objects more interesting by using combinations of randomly generated parameters. Different random generators can be specified for different parameters such as brightness and redshift. For example we can sample the brightness from a Gaussian distribution and sample the redshift from a uniform distribution.

model3 = ConstantSEDModel(
    brightness=NumpyRandomFunc("normal", loc=20.0, scale=2.0),
    redshift=NumpyRandomFunc("uniform", low=0.1, high=0.5),
    node_label="test",
)

num_samples = 10
state = model3.sample_parameters(num_samples=num_samples)
for i in range(num_samples):
    print(f"{i}: brightness={state['test']['brightness'][i]} redshift={state['test']['redshift'][i]}")

0: brightness=18.853255713574963 redshift=0.49389614504354473
1: brightness=22.319969582457624 redshift=0.47257936012114643
2: brightness=19.26396168936985 redshift=0.3247910167868396
3: brightness=23.095639400811923 redshift=0.37180460238682134
4: brightness=17.37521997677128 redshift=0.3178491137082222
5: brightness=16.911226921800136 redshift=0.17558151222426363
6: brightness=16.596162411753493 redshift=0.4503746119002576
7: brightness=18.26051019124456 redshift=0.20853632826074314
8: brightness=21.50624455807803 redshift=0.13259419859234844
9: brightness=17.839965368949084 redshift=0.4355033730627734

The sampling process creates vectors of samples for each parameter such that the i-th value of each parameter is from the same sampling run. So in the output above, sample 0 consists of all the parameter values for that sample (everything at index=0), sample 1 consists of all parameter values for that sample (everything at index=1), and so forth. This is critically important once we start dealing with parameters that are not independent. We might want to choose the host galaxy’s mass and redshift from a joint distribution. Alternatively, as we will see below, we will often want to compute one parameter as a mathematical transform of another or sample one parameter based on the value of another. In all of these cases it is important that we can access the parameters that were generated together and that these parameters stay consistent.

We can slice out a single sample using extract_single_sample() and display it.

single_sample = state.extract_single_sample(0)
print(str(single_sample))

test:
    ra: None
    dec: None
    redshift: 0.49389614504354473
    t0: None
    distance: None
    brightness: 18.853255713574963
NumpyRandomFunc:uniform_1:
    low: 0.1
    high: 0.5
    function_node_result: 0.49389614504354473
NumpyRandomFunc:normal_2:
    loc: 20.0
    scale: 2.0
    function_node_result: 18.853255713574963

You’ll notice that there are more parameters than we manually set. Some of these are for internal bookkeeping. Parameters are created automatically by the nodes if needed. In general the user should not need to worry about these extra parameters. They can access the ones of interest with the dictionary notation.

Function Nodes

Sampling functions, such as those provided by numpy, are only one type of function that we might want to use to generate parameters. We might want to sample from other functions or apply a mathematical transform to multiple input parameters to compute a new parameter. For example consider the case of computing the distmod parameter from redshift. We can do this using the information about the cosmology, such as provided by astropy’s FlatLambdaCDM class:

from astropy.cosmology import FlatLambdaCDM

cosmo_obj = FlatLambdaCDM(H0=73.0, Om0=0.3)
redshifts = np.array([0.1, 0.2, 0.3])
distmods = cosmo_obj.distmod(redshifts).value
print(distmods)

[38.22408047 39.8651567  40.86433594]

Importantly, you can use a FunctionNode that takes input parameter(s) and produces output parameter(s) during the generation. FunctionNode is a subclass of ParameterizedNode that wraps the functionality of collecting the inputs, running the computations, and storing the output. The user simply needs to provide the function node with the function it will use and the parameters. For example, LightCurveLynx has a DistModFromRedshift class that wraps the previous operation.

The below code samples a redshift uniformly from [0.1, 0.5], uses it to compute the distmod parameter, and sets that.

from lightcurvelynx.astro_utils.snia_utils import DistModFromRedshift

distmod_obj = DistModFromRedshift(
    H0=73.0, Omega_m=0.3, redshift=NumpyRandomFunc("uniform", low=0.1, high=0.5)
)

or more concretely, we can create our own FunctionNode that computes y = m * x + b.

from lightcurvelynx.base_models import FunctionNode


def _linear_eq(x, m, b):
    """Compute y = m * x + b"""
    return m * x + b


func_node = FunctionNode(
    _linear_eq,  # First parameter is the function to call.
    x=NumpyRandomFunc("uniform", low=0.0, high=10.0),
    m=5.0,
    b=-2.0,
)
print(func_node.sample_parameters(num_samples=1))

NumpyRandomFunc:uniform_1:
    low: 0.0
    high: 10.0
    function_node_result: 1.6766498498735416
FunctionNode:_linear_eq_0:
    x: 1.6766498498735416
    m: 5.0
    b: -2.0
    function_node_result: 6.383249249367708

The first parameter of the function node is the input function, such as our linear equation above. Each input into the function must be included as a named parameter, such as x, m, and b above. If any of the input parameters are missing the code will give an error.

Here we use constants for m and b so we use the same linear formulation for each sample. Only the value of x changes. However we could have also used function nodes, including sampling functions, to set m and b. In that case it is important to remember that each of our results sample will be the result of a sampling of all the variable parameters.

It would be tiresome to manually generate a FunctionNode object or class for every small mathematical function we need to use. As such LightCurveLynx also provides the BasicMathNode, which will take a string and (safely) compile the mathematical expression into a function.

from lightcurvelynx.math_nodes.basic_math_node import BasicMathNode

math_node = BasicMathNode("a + b", a=5.0, b=10.0)
print(math_node.sample_parameters(num_samples=1))

BasicMathNode:eval_0:
    a: 5.0
    b: 10.0
    function_node_result: 15.0

Linked Parameters / Nodes

Often the values of one node might depend on the values of another. A great case of this is a source/host pair where the location of the object depends on that of the host. We can access another node’s sampled parameters using a . notation: {model_object}.{parameter_name}

host = ConstantSEDModel(brightness=15.0, ra=1.0, dec=2.0, node_label="host")
source = ConstantSEDModel(brightness=10.0, ra=host.ra, dec=host.dec, node_label="source")
state = source.sample_parameters(num_samples=5)

for i in range(5):
    print(
        f"{i}: Host=({state['host']['ra'][i]}, {state['host']['dec'][i]})"
        f"Source=({state['source']['ra'][i]}, {state['source']['dec'][i]})"
    )

0: Host=(1.0, 2.0)Source=(1.0, 2.0)
1: Host=(1.0, 2.0)Source=(1.0, 2.0)
2: Host=(1.0, 2.0)Source=(1.0, 2.0)
3: Host=(1.0, 2.0)Source=(1.0, 2.0)
4: Host=(1.0, 2.0)Source=(1.0, 2.0)

We can combine the node-parameter references with functional nodes to perform actions such as sampling with noise.

Here we generate the host’s (RA, dec) from a uniform patch of the sky and then generate the source’s (RA, dec) using a Gaussian distribution centered on the host’s position. For each sample the host and source should be close, but not necessarily identical.

host = ConstantSEDModel(
    brightness=15.0,
    ra=NumpyRandomFunc("uniform", low=10.0, high=15.0),
    dec=NumpyRandomFunc("uniform", low=-10.0, high=10.0),
    node_label="host",
)

source = ConstantSEDModel(
    brightness=100.0,
    ra=NumpyRandomFunc("normal", loc=host.ra, scale=0.1),
    dec=NumpyRandomFunc("normal", loc=host.dec, scale=0.1),
    node_label="source",
)
state = source.sample_parameters(num_samples=10)

ax = plt.figure().add_subplot()
ax.plot(state["host"]["ra"], state["host"]["dec"], "b.")
ax.plot(state["source"]["ra"], state["source"]["dec"], "r.")

for i in range(5):
    print(
        f"{i}: Host=({state['host']['ra'][i]}, {state['host']['dec'][i]})    "
        f"Source=({state['source']['ra'][i]}, {state['source']['dec'][i]})"
    )

0: Host=(13.298229354895744, -1.5737770265606166)    Source=(13.19740168881626, -1.4793347389235534)
1: Host=(11.941988232359245, 8.088276450643633)    Source=(11.890974553114653, 8.270418373429516)
2: Host=(14.86831114858859, 8.921292701070854)    Source=(14.814136945626329, 9.075347205910399)
3: Host=(13.064702747958432, -6.411288200252239)    Source=(13.11080848530005, -6.4260754514450324)
4: Host=(12.026506858263097, 2.489582582457592)    Source=(12.030847401585222, 2.4513201293486113)

Again we can access all the information for a single sample. Here we see the full state tracked by the system. In addition to the host and source nodes we created, the information for the functional nodes is tracked.

single_sample = state.extract_single_sample(0)
print(str(single_sample))

NumpyRandomFunc:uniform_3:
    low: 10.0
    high: 15.0
    function_node_result: 13.298229354895744
host:
    ra: 13.298229354895744
    dec: -1.5737770265606166
    redshift: None
    t0: None
    distance: None
    brightness: 15.0
NumpyRandomFunc:uniform_4:
    low: -10.0
    high: 10.0
    function_node_result: -1.5737770265606166
NumpyRandomFunc:normal_1:
    loc: 13.298229354895744
    scale: 0.1
    function_node_result: 13.19740168881626
source:
    ra: 13.19740168881626
    dec: -1.4793347389235534
    redshift: None
    t0: None
    distance: None
    brightness: 100.0
NumpyRandomFunc:normal_5:
    loc: -1.5737770265606166
    scale: 0.1
    function_node_result: -1.4793347389235534

It is interesting to note that functional nodes themselves are parameterized nodes, allowing for more complex forms of chaining. For example we could set the low parameter from one of the NumpyRandomFuncs from another function node. This allows us to specify priors and comlex distributions.

Similarly we can now make the input parameters of one node depend on a function of parameters from other nodes. We can arbitrarily chain the computations.

# Create a host galaxy with a random brightness.
host = ConstantSEDModel(
    brightness=NumpyRandomFunc("uniform", low=1.0, high=5.0),
    node_label="host",
)

# Create the brightness of the source as a uniformly sampled foreground brightness
# added to the 80% of the host's brightness (background).
source_brightness = BasicMathNode(
    "0.8 * val1 + val2",
    val1=host.brightness,
    val2=NumpyRandomFunc("uniform", low=1.0, high=2.0),
    node_label="plus_80percent",
)

source = ConstantSEDModel(
    brightness=source_brightness,
    node_label="source",
)
state = source.sample_parameters(num_samples=10)
print(f"Host Brightness: {state['host']['brightness']}")
print(f"Source Brightness: {state['source']['brightness']}")

Host Brightness: [1.27291788 1.39326798 2.38444639 1.14544545 2.03494416 4.85321208
 2.23526331 1.41055266 1.11897484 3.94097103]
Source Brightness: [2.92036679 2.23699286 3.10984181 2.63154951 2.74103622 5.80978253
 2.94186219 2.16424567 1.99131919 4.45791253]

The state is used within the evaluate_sed() function to generate the flux densities.

time = np.arange(0, 10, 0.1)
waves = np.array([1000.0, 2000.0])
fluxes = source.evaluate_sed(time, waves, state)
print(f"Generated {fluxes.shape} fluxes (samples x times x wavelengths).")

# Plot the the flux for sample=0 and wavelength=1000.0.
plt.plot(time, fluxes[0, :, 0], "k-")
plt.xlabel("Time (days)")
plt.ylabel("Flux (nJy)")
plt.title("Flux for sample=0 and wavelength=1000.0")
plt.show()

Generated (10, 100, 2) fluxes (samples x times x wavelengths).

MultiObjectModels

We expect that many users will want to simulate fluxes produced by a combination of objects, such as a supernova and its host galaxy. LightCurveLynx provides the AdditiveMultiObjectModel for computing such combinations. The flux produced by the model is the (weighted) sum of fluxes from the individual sources.

Each underlying model in the AdditiveMultiObjectModel separately applies rest frame effects to its component flux. This allows users to model unresolved objects at different distances (with different redshifts or dust extinctions). All observer frame effects are applied to the summed fluxes for consistency.

from lightcurvelynx.models.basic_models import SinWaveModel
from lightcurvelynx.models.multi_object_model import AdditiveMultiObjectModel

# We model the host as a galaxy with a random brightness and position.
host = ConstantSEDModel(
    brightness=NumpyRandomFunc("normal", loc=10.0, scale=1.0),
    ra=NumpyRandomFunc("uniform", low=10.0, high=15.0),
    dec=NumpyRandomFunc("uniform", low=-10.0, high=10.0),
    node_label="host",
)

# We model the source as a sine wave with a given frequency and amplitude.
source = SinWaveModel(
    brightness=1.0,
    amplitude=1.0,
    frequency=0.5,
    ra=NumpyRandomFunc("normal", loc=host.ra, scale=0.1),
    dec=NumpyRandomFunc("normal", loc=host.dec, scale=0.1),
    t0=0.0,
    node_label="sin_wave_source",
)

# We combine the host and source into a multi-source model and sample it.
combined = AdditiveMultiObjectModel(
    objects=[host, source],
    node_label="combined_model",
)
state = combined.sample_parameters(num_samples=1)
print(state)

NumpyRandomFunc:uniform_2:
    low: 10.0
    high: 15.0
    function_node_result: 10.666904827138275
host:
    ra: 10.666904827138275
    dec: -7.678639148147703
    redshift: None
    t0: None
    distance: None
    brightness: 8.896906761731323
NumpyRandomFunc:uniform_3:
    low: -10.0
    high: 10.0
    function_node_result: -7.678639148147703
NumpyRandomFunc:normal_4:
    loc: 10.0
    scale: 1.0
    function_node_result: 8.896906761731323
NumpyRandomFunc:normal_6:
    loc: 10.666904827138275
    scale: 0.1
    function_node_result: 10.680212051797444
sin_wave_source:
    ra: 10.680212051797444
    dec: -7.594692205879659
    redshift: None
    t0: 0.0
    distance: None
    brightness: 1.0
    amplitude: 1.0
    frequency: 0.5
NumpyRandomFunc:normal_7:
    loc: -7.678639148147703
    scale: 0.1
    function_node_result: -7.594692205879659
combined_model:
    ra: None
    dec: None
    redshift: None
    t0: None
    distance: None

# Evaluate and generate the fluxes from the combined model for a
# single wavelength parameter sample.
time = np.arange(0, 10, 0.1)
fluxes = combined.evaluate_sed(time, np.array([1000.0]), state)

# Plot the the flux for sample=0 and wavelength=1000.0.
plt.plot(time, fluxes[:, 0], "k-")
plt.xlabel("Time (days)")
plt.ylabel("Flux (nJy)")
plt.title("Flux for combined model")
plt.show()