Simulation#

uvex_transients.simulation is where everything else in the package comes together: a transients population, sampled against a real SurveySchedule, becomes a Monte Carlo catalog of events – narrowed down to the ones that actually matter – with real synthetic photometry available for any one of them on demand. SurveySimulator does the sampling and screening; EventCatalog holds the result as a plain, portable data table; and Event reconstructs one row of that table into something you can query for a light curve or a full synthetic detection.

This page walks through that pipeline end to end. See Models and Transients for the layers underneath it, and Surveys for the schedule this all gets sampled against.

Quick Look#

The fastest way to get a feel for the pipeline is to run it. As in Surveys, we build a small synthetic schedule here (400 randomly-pointed 3x3 deg fields over 180 days) rather than downloading a real one, so the example below runs offline:

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.coordinates import EarthLocation, SkyCoord
from astropy.table import QTable
from astropy.time import Time
from regions import RectangleSkyRegion

from uvex_transients.surveys.base import SurveySchedule
from uvex_transients.simulation.core import SurveySimulator
from uvex_transients.transients.TDEs import TidalDisruptionEvent

n = 400
rng = np.random.default_rng(0)

table = QTable()
table["start_time"] = Time("2025-01-01T00:00:00") + np.sort(rng.uniform(0, 180, n)) * u.day
table["duration"] = np.full(n, 900.0) * u.s
table["observer_location"] = EarthLocation.from_geodetic(0 * u.deg, 0 * u.deg, 600 * u.km)
table["action"] = np.full(n, "observe")
table["target_coord"] = SkyCoord(
    rng.uniform(0, 360, n) * u.deg,
    np.degrees(np.arcsin(rng.uniform(-1, 1, n))) * u.deg,
)
table["roll"] = np.zeros(n) * u.deg
table["field_id"] = np.arrange(n)
table["block_id"] = np.zeros(n, dtype=int)

fov = RectangleSkyRegion(center=SkyCoord(0 * u.deg, 0 * u.deg), width=3 * u.deg, height=3 * u.deg)
schedule = SurveySchedule(table, fov)

tde = TidalDisruptionEvent()
simulator = SurveySimulator(schedule, transients={"tde": tde}, simulation_seed=0)

catalog = simulator.generate_events(time_bins=6, nside=32)

fig = plt.figure(figsize=(7, 4))
ax = fig.add_subplot(111, projection="aitoff")
ax.grid(True)
ra = catalog.coord.ra.wrap_at(180 * u.deg).radian
ax.scatter(ra, catalog.coord.dec.radian, s=4, color="C0")
ax.set_title(f"{len(catalog)} sampled TDEs across the example schedule")

(Source code)

That’s the whole shape of it: pair a transient population with a schedule, hand both to a generate_events() call, and get back every TDE that could plausibly have exploded somewhere the schedule looked. The rest of this page covers the simulator itself, the two screening steps that narrow that population down to detections, and how to turn any one surviving row back into a real light curve.

Hint

time_bins above is a plain integer – the number of evenly-spaced bins to divide the schedule’s own span into. It also accepts an explicit Time array of bin edges, for windows that don’t line up evenly with the schedule (e.g. one bin per lunation).

The Simulator Object#

SurveySimulator is a thin coordinator: a SurveySchedule, a {name: transient} dict of ExtragalacticTransient instances to sample, and a root seed for reproducibility, all supplied at construction:

from uvex_transients.simulation.core import SurveySimulator
from uvex_transients.transients.TDEs import TidalDisruptionEvent
from uvex_transients.transients.kilonovae import Kilonova

simulator = SurveySimulator(
    schedule,
    transients={"tde": TidalDisruptionEvent(), "kne": Kilonova()},
    simulation_seed=42,
)

simulator.survey_schedule       # the SurveySchedule passed in
simulator.transient_collection  # the {name: transient} dict
simulator.simulation_seed       # 42

Important

Every transient in transients needs its own distinct name – this is the key that ties each sampled event back to the transient instance that produced it (the catalog’s transient_type column) and, later, that reconstructs an Event from a catalog row via get_events(). See Transients for configuring a transient’s priors, duration_limit, or redshift limit before handing it to the simulator.

Simulating several transient types together, as above, is no different from simulating one – every step below (generate_events, both filtering methods, and EventCatalog.get_events) loops over every transient_type present in the catalog, keyed against this same transient_collection dict.

Generating the Event Catalog#

generate_events() is the main event: a Monte Carlo realization of every registered transient type, sampled only where and when the schedule could plausibly have caught it.

Windowed Sampling#

Sampling a transient’s full redshift- and sky-dependent volumetric rate over the entire sky and survey duration, then throwing away everything the schedule never observed, would waste almost all of that effort – a real survey only ever covers a small fraction of the sky at any one time. To avoid that, generate_events() divides the schedule into time_bins windows and, for each window and each transient type, first asks get_observed_healpix_ids() which HEALPix pixels the schedule touches at all between the start of the window and its end plus that transient’s own duration_limit – i.e. late enough that a transient exploding right at the end of the window could still be caught while it’s active. Events are then sampled only within those pixels, via sample_events_on_healpix_grid(), with explosion times drawn uniformly within the window itself (never the padded tail), so no event is ever double-counted across adjacent bins. See Transients for how that per-pixel sampling itself works.

catalog = simulator.generate_events(time_bins=6, nside=32, order="nested")

nside/order set the HEALPix resolution both the coverage lookup and the sampling grid use – the same tradeoff as everywhere else HEALPix shows up in this package: finer pixels track the schedule’s real footprint more closely, at the cost of more pixels to sample per bin.

Tip

For a transient type with a very high intrinsic rate, downsample=k draws a random 1/k subset of each per-bin, per-type table (without replacement, seeded off simulation_seed) instead of sampling the full population. Multiply any downstream count by k to get back an estimate of the true yield; see the Simulating Transients example for this in practice. downsample can instead be a {transient key: k} mapping to downsample types individually – a type left out of the mapping isn’t downsampled at all. Either form is stashed on the returned EventCatalog as downsample, purely as provenance (nothing rescales counts back up automatically), and every cut carries it through to its own output catalog unchanged.

Two columns are computed once here, rather than being left for every later step to re-derive: each event’s luminosity_distance (interpolated off that transient type’s own cached luminosity_distance_grid, not a fresh cosmology call per event) and its ebv (one vectorized Milky Way foreground dust-map query, via uvex_transients.dust, over every sampled position at once).

The Event Catalog#

generate_events() returns an EventCatalog: a plain data table with no live reference back to the schedule or the transient instances it was generated from (both are supplied again, explicitly, wherever they’re needed – see Reconstructing and Simulating Events below). That makes it trivially picklable and safe to round-trip to disk:

len(catalog)          # number of sampled events
catalog.table          # the underlying astropy.table.QTable

catalog.to_disk("tde_catalog.ecsv", overwrite=True)
reloaded = EventCatalog.from_disk("tde_catalog.ecsv")

to_disk() writes the table as ECSV with nside/order/time_bins/seed/downsample stashed in the file’s header, so from_disk() can reconstruct a complete EventCatalog from the one file alone. (A file written by an older version of the package, with no downsample in its header, still reads back fine – it just defaults to None.)

Every column of catalog.table is also available as a convenience property, returning a plain array (or Quantity/Time/ SkyCoord, as appropriate) rather than a raw table column:

Column / property	Type	Meaning
`event_id`	int	Unique id, assigned once across the whole catalog (not renumbered by filtering).
`transient_type`	str	Which entry of `transient_collection` this event came from.
`time_bin`	int	Index into `time_bins` of the window the event exploded within.
`healpix_id`, `healpix_dx`, `healpix_dy`	int, float, float	Sampling pixel (at this catalog’s own `nside`/`order`) and the event’s sub-pixel offset within it.
`coord`	`SkyCoord`	Sky position.
`redshift`	float	Cosmological redshift.
`luminosity_distance`	`Quantity`	Cached at generation time; see Simulation above.
`ebv`	float	Milky Way foreground E(B-V) at `coord`; cached at generation time.
`t_explosion`	`Time`	Time of explosion.
`parameter_seed`	int	Regenerates this event’s physical SED parameters on demand – see Reconstructing and Simulating Events below.

Note

An event’s physical SED parameters (amplitude, rise time, temperature, …) are deliberately not stored as columns – only the parameter_seed that regenerates them. Storing one column per parameter would mean a different schema per transient type; regenerating them lazily from a stored seed keeps EventCatalog – and every filtering method below – agnostic to which SED any given row actually uses.

Filtering the Event Catalog#

A freshly-sampled catalog is dominated by events far too faint to ever matter – most of a transient’s redshift-limited volume is, by construction, near the limit where it’s essentially undetectable. Two progressively more expensive passes narrow it down to the events actually worth keeping, both taking an EventCatalog and an Mission (for its Detector’s bandpasses) and returning a new EventCatalog over the surviving rows – nside/order/time_bins/seed/downsample unchanged, and original event_id values preserved rather than renumbered:

Limiting Magnitude

Method: filter_by_limiting_magnitude()

Description: Schedule-independent – it never checks whether the schedule actually pointed anywhere near an event. For each transient type, every event’s own physical parameters are regenerated from its stored parameter_seed and evaluated over a shared linspace(0, duration_limit, n_phase) phase grid, the same grid for every event of that type; an event survives if its brightest requested band clears mag_limit at n_visits or more of those phase samples.

mag_filtered = simulator.filter_by_limiting_magnitude(
    catalog, mission, mag_limit=25.0,
)

Uses: A cheap first pass: “could this event, at its absolute brightest, ever be seen at all?” chunk_size bounds peak memory (evaluation broadcasts every event and phase sample into one dense array at once), and raising n_visits above its default of 1 discards events that only momentarily clear the limit.

Signal-to-Noise Ratio

Method: filter_by_snr()

Description: The real question, schedule and all: over every observation the schedule actually made of an event’s position while it was active, is it ever detected above snr_threshold? Every surviving event’s observations are gathered with one batched get_observation_indices_of() call per chunk (see Surveys), then evaluated together with get_snr().

detected = simulator.filter_by_snr(
    mag_filtered, mission, snr_threshold=5.0,
)

Uses: The actual detected population. Run it against the output of filter_by_limiting_magnitude, not the raw catalog – it’s the more expensive of the two (every remaining event costs at least one schedule query), so there’s no reason to pay that cost on events the cheap pass would have rejected anyway.

from m4opt.missions import uvex

mission = uvex
mag_filtered = simulator.filter_by_limiting_magnitude(catalog, mission, mag_limit=25.0)
detected = simulator.filter_by_snr(mag_filtered, mission, snr_threshold=5.0)

stages = ["Sampled", "Mag < 25", "SNR > 5"]
counts = [len(catalog), len(mag_filtered), len(detected)]

fig, ax = plt.subplots()
ax.bar(stages, counts, color=["#888888", "#4C72B0", "#55A868"])
for i, count in enumerate(counts):
    ax.text(i, count, f"{count:,}", ha="center", va="bottom")
ax.set_ylabel("Number of TDEs")
ax.set_title("TDE detection funnel")

(Source code)

Exposure and Yield#

A raw or filtered EventCatalog count is a realization, not an estimate – to turn one into a formal expected-detection number with confidence bounds (Yield Statistics derives every estimator below), you also need to know how much of the sky the survey actually covered. compute_effective_exposure() answers that, independent of any Monte Carlo draw: it reruns exactly the same per-bin, per-type footprint query generate_events restricts its own sampling to, but reduces it to a solid angle instead of a drawn population:

exposure = simulator.compute_effective_exposure(time_bins=6, nside=32)

exposure.total_effective_exposure   # {transient type: total solid-angle*time exposure}
exposure.total_expected_events      # {transient type: mu_0, the footprint-aware expected count}
exposure.coverage_fraction          # {transient type: fraction of the full 4*pi sky-time swept}

The returned ExposureCatalog has one row per (transient type, time bin) – the same time_bins/nside/order as generate_events should always be passed here too, so both describe the same footprint query. Each row’s expected_events is effective_exposure * transient.integrated_rate – Yield Statistics’s \(\mu_0\) for that bin – so summing it over every bin (total_expected_events) gives the intrinsic expected count actually reachable by this schedule’s footprint, not compute_all_sky_yield()’s idealized full-sky number. get_exposure_between()/ get_expected_events_between() and rebin() let you query or re-tile that same tabulated exposure over an arbitrary sub-window or a different time binning without re-querying the schedule.

Combine a raw (feasible) catalog, a detected (post-cut) catalog, and its exposure into one per-transient-type summary with compute_yield_summary():

yields = catalog.compute_yield_summary(detected, exposure, {"tde": tde}, confidence=0.9)
yields.table["transient_type", "uvex_intrinsic_events", "detection_probability", "expected_detections"]

The returned YieldTable carries, per transient type, the rate (integrated_rate/all_sky_rate), the footprint-aware intrinsic rate/count (uvex_intrinsic_rate/uvex_intrinsic_events, i.e. \(\mu_0\)), the Monte Carlo detection efficiency (detection_probability, \(\hat\epsilon=k/n\)), and the final yield estimate (expected_detections, \(\hat\lambda=\mu_0\hat\epsilon\)) – each rate-derived column with its own RATE_CI-propagated bounds, and detection_probability/expected_detections with two separate uncertainty sources (Clopper-Pearson binomial and rate-normalization), kept apart as ..._binom_lower/_upper and ..._rate_lower/_upper columns rather than combined into one. to_ascii() writes a human-readable summary table; to_latex() renders both uncertainty sources as stacked LaTeX superscripts for a paper table.

Hint

For a finer-grained question than “was this event detected at all” – “how many separate epochs was it detected in” – run synthetic photometry over a whole catalog with compute_photometry_catalog() (wrapping simulate_photometry(), below, as a PhotometryCatalog), then call its own compute_detection_count_table(). It generalizes compute_yield_summary’s “detected at all” (\(N_{\rm det}\geq 1\)) to “detected in at least \(k\) epochs” for every \(k\) at once, with the same two-source uncertainty treatment – exactly what the detection-counts CLI command (Command-Line Interface) automates.

Reconstructing and Simulating Events#

Every screening step above works with cheap, batched approximations – a shared phase grid, a single flux point at each band’s pivot wavelength. Getting a real, per-observation synthetic light curve for one particular event means reconstructing it as a full Event, via get_events(). Because an EventCatalog holds no live reference to the schedule or the transients it came from (see above), both are supplied again here – typically the same ones the catalog was generated with:

event = detected.get_events(19, {"tde": tde}, schedule)
print(event)
... <Event id=19 type='tde' z=0.7097 n_observations=2>

Building an Event runs exactly one query against the schedule (get_observations_of()) to find which scheduled observations actually covered it while active – available as observations (n_observations) – alongside its coord, redshift, luminosity_distance, ebv, and t_explosion. No photometry is done yet at this point; sample_parameters() regenerates its physical SED parameters (deterministically, from the same stored parameter_seed) as a pure, idempotent lookup any of the methods below can call as many times as needed.

Theoretical Curves#

mag(), flux(), and luminosity() give the noiseless truth this event’s real photometry (below) scatters around – the first two folding in this event’s own redshift, distance, and foreground dust; the last a rest-frame, distance-independent bolometric quantity:

t = np.linspace(0, tde.duration_limit.to_value(u.day), 300) * u.day

event.mag(t, mission, band="NUV")     # apparent AB magnitude
event.flux(t, mission, band="NUV")    # observed flux density, at the band's pivot wavelength
event.luminosity(t)                   # rest-frame bolometric L_bol(t)

Synthetic Photometry#

simulate_photometry() is the expensive step every earlier screening pass was designed to defer: real, per-observation, per-band synthetic photometry against every observation in observations, batched into one get_snr() call per band regardless of how many observations there are:

event = detected.get_events(19, {"tde": tde}, schedule)

phot = event.simulate_photometry(mission)
t_since_explosion = (phot["obs_time"] - event.t_explosion).to(u.day)
t_theory = np.linspace(0, tde.duration_limit.to_value(u.day), 300) * u.day

fig, ax = plt.subplots(figsize=(7, 4))
for band, color in {"FUV": "#4C72B0", "NUV": "#DD8452"}.items():
    ax.plot(t_theory.value, event.mag(t_theory, mission, band=band).value, color=color, lw=1.5, alpha=0.6)

    in_band = np.isfinite(phot["ab_mag"]) & (phot["band"] == band)
    if np.any(in_band):
        ax.errorbar(
            t_since_explosion[in_band].value,
            phot["ab_mag"][in_band],
            yerr=phot["mag_err"][in_band],
            marker="s", mfc=color, mec="k", ecolor=color, linestyle="none", label=band,
        )

ax.invert_yaxis()
ax.set_xlabel("Days since explosion")
ax.set_ylabel("AB magnitude")
ax.set_title(f"Event {event.event_id} (z={event.redshift:.2f}, {event.n_observations} observations)")
ax.legend()

(Source code)

This example schedule only ever visits a given field once, so most of its events – like this one – end up with just one or two associated observations; see the Simulating Transients for the same workflow against a real, densely-cadenced UVEX schedule.

Each returned row is a Gaussian realization of the true flux at that observation’s implied SNR, not the ground truth itself – flux/ab_mag and their symmetric-in-flux n_sigma bounds (flux_upper/flux_lower, transformed separately to mag_upper/mag_lower rather than through a single linearized mag_err) are nan wherever the noisy draw itself isn’t securely above zero flux – the correct behavior at low SNR, not a bug. The whole event replays identically given the same parameter_seed, since every parameter draw and every band’s noise realization comes from one RNG seeded from it.

See the Simulating Transients for a full worked example of this entire pipeline against a real UVEX schedule, and uvex_transients.simulation in the API reference for exhaustive method-by-method detail.