sys_mapping.simulation

Systematic contamination injection, FITS I/O, and w(θ) recovery.

Provides the end-to-end infrastructure for simulation-based validation:

Load catalogs — load_uchuu_mock() reads the Uchuu lightcone FITS files; GLASS catalogs come from sys_mapping.glass_mocks.
Load templates — load_systematic_maps() reads LSDR10 HEALPix systematic maps.
Build contamination grid — make_contamination_grid() returns 9 ContaminationConfig objects spanning 3 amplitude levels × 3 contamination scenarios.
Inject contamination — inject_systematics() computes per-galaxy weights \(\texttt{WEIGHT\_CONT}(p) = (1+\delta_{\rm cont}(p))/(1+\delta_g(p))\).
Save / load FITS — save_simulation_catalog() and load_simulation_catalog().
Full recovery pipeline — run_wtheta_recovery() measures truth, contaminated, and recovered \(w(\theta)\) for OLS, ISD-1, and ElasticNet.

Systematic contamination injection, FITS I/O, and w(θ) recovery.

This module provides the end-to-end infrastructure for simulation-based validation of the decontamination pipeline:

Load reference catalogs (Uchuu lightcone or GLASS mock).
Load LSDR10 systematic maps.
Inject contamination into catalog positions using the pixel-level forward model from sys_mapping.contamination.
Save contaminated catalogs to FITS.
Run all fast decontamination methods and compare recovered w(θ) to truth.

Contamination injection strategy

For a catalog of N galaxies, the true overdensity at each pixel is estimated from the catalog itself (galaxy counts / random counts − 1). The forward contamination model (Eq. 11-13 of Berlfein et al. 2024) is applied to yield a contaminated overdensity field. The per-galaxy contamination weight is then:

WEIGHT_CONT(p) = (1 + δ_cont(p)) / (1 + δ_g(p))

When computing w(θ) with TreeCorr, passing WEIGHT_CONT as galaxy weights exactly reproduces the effect that the contamination model would have on the angular correlation function, without physically removing or adding galaxies. The “truth” w(θ) is computed with uniform weights (no contamination).

References

Berlfein et al. 2024, MNRAS 531, 4954. arXiv:2401.12293

class sys_mapping.simulation.ContaminationConfig(level, scenario, a_true, b_true)[source]

Bases: object

Parameters for a single contamination scenario.

Parameters:

level (str)
scenario (str)
a_true (ndarray)
b_true (ndarray)

level

Human-readable amplitude label: "low", "medium", or "high".

Type:: str

scenario

Contamination model: "additive", "multiplicative", or "combined".

Type:: str

a_true

Additive amplitudes (shape (n_sys,)). Zero for multiplicative.

Type:: numpy.ndarray

b_true

Multiplicative amplitudes (shape (n_sys,)). Zero for additive.

Type:: numpy.ndarray

level: str

scenario: str

a_true: ndarray

b_true: ndarray

property n_sys: int

to_header_dict()[source]

Flat dict suitable for storage in a FITS header.

Return type:: dict

sys_mapping.simulation.make_contamination_grid(n_sys, seed=0)[source]

Build a 3×3 grid of contamination configurations.

Returns 9 ContaminationConfig objects: one for each combination of three amplitude levels (low/medium/high) × three scenarios (additive/multiplicative/combined).

The signs of the amplitudes are drawn once from the given seed and shared across levels so that the impact of amplitude magnitude can be cleanly isolated.

Parameters:

n_sys (int) – Number of systematic templates.
seed (int) – Random seed for the sign draws.

Return type:

List of 9 ContaminationConfig objects.

Examples

>>> from sys_mapping.simulation import make_contamination_grid
>>> cfgs = make_contamination_grid(n_sys=3)
>>> len(cfgs)
9
>>> {c.level for c in cfgs} == {'low', 'medium', 'high'}
True
>>> {c.scenario for c in cfgs} == {'additive', 'multiplicative', 'combined'}
True

sys_mapping.simulation.load_uchuu_mock(data_fits, rand_fits)[source]

Load Uchuu lightcone FITS catalog.

Parameters:

data_fits (str | Path) – Path to the *_DATA.fits file. Expected columns: RA, DEC, redshift_S.
rand_fits (str | Path) – Path to the *_RAND.fits file. Expected columns: RA, DEC.

Returns:

dict with keys ra, dec, z (galaxies) and ra_rand,
dec_rand (randoms).

Return type:

dict

Examples

>>> import os
>>> p = os.path.expanduser(
...     '~/data/Uchuu/FullSky/mock_catalogues/'
...     'MOCK_VLIM_ANY_10.65_Mstar_12.0_0.05_z_0.26_N_0923373/'
...     'MOCK_VLIM_ANY_10.65_Mstar_12.0_0.05_z_0.26_N_0923373_DATA.fits')
>>>
>>> cat = load_uchuu_mock(p, p.replace('_DATA', '_RAND'))

sys_mapping.simulation.load_systematic_maps(syst_dir, nside, map_names=None)[source]

Load LSDR10 systematic maps at a given HEALPix resolution.

Parameters:

syst_dir (str | Path) – Directory containing NSIDE subdirectories, e.g. ~/data/legacysurvey/dr10/systematics/. The function looks for files in {syst_dir}/{nside:04d}/ (NSIDE in 4-digit zero-padded subdirectory).
nside (int) – Target HEALPix resolution.
map_names (list[str] | None) – Which maps to load. Defaults to DEFAULT_MAP_NAMES (5 maps).

Returns:

templates ((n_sys, 12*nside²) array, zero-mean and unit-std over valid pixels.)
names (list of map name strings.)
footprint ((12*nside²,) bool array — True where *all loaded maps are valid.*) – Use this to restrict mock catalogs to the survey footprint before injecting systematics or measuring w(θ).

Return type:

tuple[ndarray, list[str], ndarray]

Examples

>>>
>>> t, names, footprint = load_systematic_maps('~/data/legacysurvey/dr10/systematics/', 64)
>>> t.shape
(5, 49152)
>>> footprint.sum()  # number of pixels inside the LS10 footprint
...

sys_mapping.simulation.apply_footprint_mask(catalog, footprint, nside)[source]

Filter a catalog dict to positions inside the survey footprint.

Keeps only galaxies and randoms whose HEALPix pixel (at nside) is flagged as valid in footprint. All other keys in catalog are passed through unchanged.

Parameters:

catalog (dict) – Dict with at least ra, dec, ra_rand, dec_rand keys (and optionally z). As returned by generate_glass_fullsky_mock().
footprint (ndarray) – (12*nside²,) bool array — True = inside survey footprint. As returned by load_systematic_maps().
nside (int) – HEALPix NSIDE matching footprint.

Returns:

Filtered catalog dict. n_total (if present) is updated to the number
of galaxies remaining after filtering.

Return type:

dict

Notes

This must be called before inject_systematics() and run_wtheta_recovery() so that w(θ) is measured only within the footprint where systematic templates are defined.

sys_mapping.simulation.inject_systematics(ra, dec, ra_rand, dec_rand, templates, a, b, nside)[source]

Compute per-galaxy contamination weights.

The function estimates the pixel-level galaxy overdensity from the input catalog, applies the forward contamination model, and returns per-galaxy weights that reproduce the contaminated angular correlation function when passed to TreeCorr.

Parameters:

ra (ndarray) – Galaxy positions in degrees.
dec (ndarray) – Galaxy positions in degrees.
ra_rand (ndarray) – Random catalog positions in degrees.
dec_rand (ndarray) – Random catalog positions in degrees.
templates (ndarray) – Systematic template maps (shape (n_sys, n_pix)).
a (ndarray) – Additive contamination amplitudes (shape (n_sys,)).
b (ndarray) – Multiplicative contamination amplitudes (shape (n_sys,)).
nside (int) – HEALPix NSIDE used for pixelisation.

Returns:

WEIGHT_CONT ((n_gal,) per-galaxy weights. Equal to 1 when a=b=0.)
Clipped to a minimum of 0.01 to prevent unphysical negative weights.

Return type:

ndarray

Notes

The weight is defined as:

WEIGHT_CONT(p) = (1 + δ_cont(p)) / (1 + δ_g(p))

where δ_cont is the contaminated overdensity and δ_g is the observed overdensity from the input catalog. Pixels where 1 + δ_g(p) ≤ 0 are assigned weight 1.

Examples

>>> import numpy as np
>>> from sys_mapping.simulation import inject_systematics
>>> rng = np.random.default_rng(0)
>>> n_gal, n_rand = 5000, 50000
>>> ra = rng.uniform(0, 360, n_gal)
>>> dec = rng.uniform(-90, 90, n_gal)
>>> ra_r = rng.uniform(0, 360, n_rand)
>>> dec_r = rng.uniform(-90, 90, n_rand)
>>> templates = rng.standard_normal((2, 12 * 32**2))
>>> a = np.array([0.0, 0.0])
>>> b = np.array([0.0, 0.0])
>>> w = inject_systematics(ra, dec, ra_r, dec_r, templates, a, b, nside=32)
>>> w.shape == (n_gal,)
True
>>> np.allclose(w, 1.0, atol=0.15)   # near unity for zero contamination
True

sys_mapping.simulation.save_simulation_catalog(ra, dec, z, weight_cont, config, output_path)[source]

Write a contaminated simulation catalog to a FITS BinTableHDU.

Parameters:

ra (ndarray) – Galaxy positions in degrees.
dec (ndarray) – Galaxy positions in degrees.
z (ndarray) – Galaxy redshifts.
weight_cont (ndarray) – Per-galaxy contamination weights (from inject_systematics).
config (ContaminationConfig) – Contamination configuration stored in the FITS header.
output_path (str | Path) – Destination FITS file path. Parent directory is created if needed.

Return type:

None

sys_mapping.simulation.load_simulation_catalog(fits_path)[source]

Load a contaminated simulation catalog from FITS.

Return type:: ra, dec, z, weight_cont, config
Parameters:: fits_path (str | Path)

sys_mapping.simulation.run_wtheta_recovery(catalog, templates, template_names, config, nside=64, methods=None, min_sep=0.1, max_sep=10.0, nbins=10)[source]

Run full decontamination pipeline and return w(θ) at each stage.

Given a clean catalog and contamination parameters, this function:

Computes the reference (truth) w(θ) from the clean catalog.
Injects contamination to build per-galaxy WEIGHT_CONT.
Computes the contaminated w(θ) using those weights.
For each decontamination method, fits the contamination model and applies per-galaxy correction weights; computes the recovered w(θ).

Parameters:

catalog (dict) – Dict with keys ra, dec, z, ra_rand, dec_rand.
templates (ndarray) – Full-sky systematic maps (shape (n_sys, n_pix)).
template_names (list[str]) – Names of the template maps (for metadata in output).
config (ContaminationConfig) – Contamination configuration to inject.
nside (int) – HEALPix resolution for overdensity computation.
methods (list[str] | None) – Decontamination methods to run. Default: OLS, ISD-1, ElasticNet.
min_sep (float) – Angular separation range in degrees.
max_sep (float) – Angular separation range in degrees.
nbins (int) – Number of angular bins.

Returns:

dict with keys
- ``theta`` (bin centres in degrees)
- ``w_true`` (w(θ) of clean catalog (truth))
- ``w_contaminated`` (w(θ) with WEIGHT_CONT applied)
- ``w_recovered_{method}`` (w(θ) after decontamination for each method)
- ``params_{method}`` (dict of recovered a_hat, b_hat)
- ``config`` (the ContaminationConfig used)
- ``n_gal``, ``n_rand`` (catalog sizes)

Return type:

dict