sys_mapping.diagnostics

Post-correction diagnostics: null tests, SNR-based template ranking, and footprint masking sensitivity analysis.

Inputs: Per-pixel weight array weights, template maps delta_t, observed overdensity delta_g_obs.

Outputs: Correlation arrays, p-values, SNR arrays, or dictionaries of masking-level results.

• null_test_cross_correlations() — Pearson \(r(w, t_i)\) with permutation p-values; should be near zero after correction.

• snr_template_ranking() — three SNR estimators ("template", "data", "peak") to rank which templates carry the most contaminating power.

• footprint_mask_diagnostics() — stability of fitted amplitudes under varying mask thresholds.

Key papers: Ross et al. 2011; Tanidis et al. 2026; Weaverdyck & Huterer 2021; Rodríguez-Monroy et al. 2025 — see also Methods Reference.

Systematic mapping diagnostics: null tests, SNR ranking, footprint masking.

Implements:

• Null test cross-correlations (Ross+2011)

• SNR-based template ranking with three estimators (Tanidis+2026)

• Footprint masking sensitivity diagnostics (Rodríguez-Monroy+2025)

sys_mapping.diagnostics.null_test_cross_correlations(weights, delta_t, n_bootstrap=100, seed=0)

Test for residual systematic contamination via weight × template correlations.

After decontamination, the per-pixel systematic weights should be statistically independent of all template maps. A significant Pearson correlation between weights and delta_t[i] indicates residual contamination from template i.

Parameters:

• weights (ndarray) – Per-pixel systematic weights (shape (n_pix,)).

• delta_t (ndarray) – Template maps (shape (n_sys, n_pix)).

• n_bootstrap (int) – Number of permutation resamples for computing p-values.

• seed (int) – Random seed for reproducibility.

Returns:

Dictionary with:

• "correlations": Pearson \(r(w, t_i)\) for each template (shape (n_sys,)).

• "p_values": two-sided permutation p-value for each template (shape (n_sys,)). Small values indicate significant residual contamination.

Return type:

result

Notes

The Pearson correlation is:

\[r_i = \frac{\sum_p (w_p - \bar w)(t_{i,p} - \bar t_i)} {\sqrt{\sum_p (w_p-\bar w)^2 \sum_p (t_{i,p}-\bar t_i)^2}}\]

A well-corrected field satisfies \(|r_i| \approx 0\) for all templates. Permutation p-values are computed by shuffling weights and recomputing the correlation n_bootstrap times.

Examples

>>> import numpy as np
>>> from sys_mapping import null_test_cross_correlations
>>> rng = np.random.default_rng(0)
>>> n_pix, n_sys = 500, 3
>>> weights = rng.standard_normal(n_pix) + 1.0
>>> delta_t = rng.standard_normal((n_sys, n_pix))
>>> result = null_test_cross_correlations(weights, delta_t, n_bootstrap=50)
>>> result["correlations"].shape
(3,)
>>> result["p_values"].shape
(3,)
>>> np.all(result["p_values"] >= 0) and np.all(result["p_values"] <= 1)
True

References

Ross et al. 2011, MNRAS 417, 1350.

sys_mapping.diagnostics.snr_template_ranking(delta_g_obs, delta_t, *, method='template', n_bins=10, poly_order=1, fracdet=None)

Rank systematic templates by signal-to-noise ratio of their contamination.

Four SNR definitions:

• "template": \(|\hat\alpha_i| / \sigma_{\hat\alpha_i}\) from OLS.

• "data": \(|\mathrm{Corr}(\delta_g, t_i)|\), the absolute cross-correlation between the galaxy field and each template.

• "peak": peak pseudo-\(C_\ell\) cross-spectrum between \(\delta_g\) and \(t_i\), divided by the noise level.

• "isd": \(\Delta\chi^2 = \chi^2_{\rm null} - \chi^2_{\rm model}\) from the ISD 1D binned relation (Rodríguez-Monroy et al. 2025, Eqs. 4–5). For each template, pixels are binned by template value; the fracdet-weighted mean galaxy density per bin is compared to a polynomial fit. Large \(\Delta\chi^2\) indicates significant systematic contamination.

Parameters:

• delta_g_obs (ndarray) – Observed galaxy overdensity (shape (n_pix,)).

• delta_t (ndarray) – Template maps (shape (n_sys, n_pix)).

• method (str) – One of "template", "data", "peak", "isd".

• n_bins (int) – Number of equal-width bins for method="isd". Ignored otherwise.

• poly_order (int) – Polynomial degree for the 1D fit in method="isd" (1 = linear, 3 = cubic as used in DES Y6). Ignored otherwise.

• fracdet (ndarray | None) – Per-pixel fractional coverage weights (shape (n_pix,)). Only used by method="isd"; if None, uniform weights are assumed.

Returns:

SNR value for each template (shape (n_sys,)). Higher values indicate a more contaminating template.

Return type:

snr

Notes

Templates can be sorted by np.argsort(snr)[::-1] to obtain a ranking from most to least contaminating. The "peak" method requires healpy (already a package dependency) and assumes the overdensity map covers the full sky (or has zeros outside the mask).

Examples

>>> import numpy as np
>>> from sys_mapping import snr_template_ranking
>>> rng = np.random.default_rng(0)
>>> n_pix, n_sys = 3000, 4
>>> delta_t = rng.standard_normal((n_sys, n_pix))
>>> # Inject strong contamination from template 2
>>> delta_g = 0.5 * delta_t[2] + rng.standard_normal(n_pix) * 0.2
>>> snr = snr_template_ranking(delta_g, delta_t, method="data")
>>> snr.shape
(4,)
>>> np.argmax(snr) == 2
True

References

Tanidis et al. 2026, MNRAS 547. Rodríguez-Monroy et al. 2025, arXiv:2509.07943, Sec. IV.A.1, Eqs. 4–5.

sys_mapping.diagnostics.footprint_mask_diagnostics(delta_g_obs, delta_t, mask_fractions, good_pixels)

Quantify sensitivity of systematic parameters to the masking threshold.

For each masking fraction in mask_fractions, the most poorly connected pixels (those with the lowest overdensity variance from the template model) are excluded, and OLS contamination amplitudes are re-estimated. Large variation of the amplitudes with masking threshold indicates sensitivity to the survey boundary or to systematic outliers.

This is the footprint masking sensitivity analysis of Rodríguez-Monroy et al. 2025, Sec. 4.

Parameters:

• delta_g_obs (ndarray) – Observed galaxy overdensity at unmasked pixels (shape (n_pix,)).

• delta_t (ndarray) – Template maps at unmasked pixels (shape (n_sys, n_pix)).

• mask_fractions (ndarray) – Array of additional pixel fractions to mask on top of the baseline, in increasing order (e.g. [0.05, 0.10, 0.15, 0.20]).

• good_pixels (ndarray) – Boolean mask over the full footprint (shape (n_full_pix,)). Used only for computing a pixel-quality ranking (template variance at each pixel).

Returns:

Dictionary with:

• "alpha_hat": OLS amplitudes per masking level (shape (n_thresholds, n_sys)).

• "scatter": standard deviation of alpha_hat across masking levels (shape (n_sys,)). Near-zero scatter means the amplitude estimates are robust to masking.

Return type:

result

Notes

Pixels are ranked by the RMS template value \(\sigma_p = \sqrt{\sum_i t_{i,p}^2 / n_{\rm sys}}\). Pixels with the lowest \(\sigma_p\) (fewest systematic features) are masked first, so that the remaining pixels have progressively higher systematic leverage.

Examples

>>> import numpy as np
>>> from sys_mapping import footprint_mask_diagnostics
>>> rng = np.random.default_rng(0)
>>> n_pix, n_sys = 4000, 3
>>> delta_t = rng.standard_normal((n_sys, n_pix))
>>> delta_g = 0.1 * delta_t[0] + rng.standard_normal(n_pix) * 0.5
>>> good = np.ones(n_pix, dtype=bool)
>>> fracs = np.array([0.0, 0.05, 0.10, 0.15])
>>> result = footprint_mask_diagnostics(delta_g, delta_t, fracs, good)
>>> result["alpha_hat"].shape
(4, 3)
>>> result["scatter"].shape
(3,)

References

Rodríguez-Monroy et al. 2025, arXiv:2509.07943.

sys_mapping.diagnostics.isd_template_significance(delta_g_obs, delta_t, good_pixels, nside, n_total, z_edges, nz, *, n_total_footprint=None, n_mocks=100, n_bins=10, poly_order=1, fracdet=None, seed=0, rand_factor=10, n_jobs=1)

ISD Δχ² significance: compare data against systematic-free GLASS mocks.

For each template map, computes the ISD contamination metric \(\Delta\chi^2 = \chi^2_{\rm null} - \chi^2_{\rm model}\) on the data and on n_mocks systematic-free GLASS mocks generated on the same footprint and redshift range. The fraction of mocks that exceed the data value defines the template’s p-value.

Parameters:

• delta_g_obs (ndarray) – Observed galaxy overdensity at footprint pixels (shape (n_good_pix,)).

• delta_t (ndarray) – Template maps at footprint pixels (shape (n_sys, n_good_pix)).

• good_pixels (ndarray) – Boolean footprint mask over the full HEALPix sky (shape (12 * nside**2,)).

• nside (int) – HEALPix resolution of the template and galaxy maps.

• n_total (int) – Target galaxy count per GLASS mock (full-sky). Ignored when n_total_footprint is provided.

• n_total_footprint (int | None) – Number of galaxies in the survey footprint (i.e. len(ra_gal) from the real data). When provided, n_total is computed automatically as n_total_footprint × n_full_pix / n_good_pix so that after trimming the mock to good_pixels, the footprint surface density matches the data. Preferred over passing a raw n_total.

• z_edges (ndarray) – Redshift bin edges for the n(z) model (length n_bins_z + 1).

• nz (ndarray) – Galaxy counts per redshift bin (length n_bins_z).

• n_mocks (int) – Number of systematic-free mocks to generate.

• n_bins (int) – Number of template-value bins for the ISD 1D relation.

• poly_order (int) – Polynomial degree for the 1D fit (1 = linear, 3 = cubic).

• fracdet (ndarray | None) – Per-pixel fractional coverage weights (shape (n_good_pix,)). If None, uniform weights are used.

• seed (int) – Base random seed; mock k uses seed + k.

• rand_factor (int) – Ratio of random to galaxy counts in each GLASS mock. Reducing this (e.g. to 2) cuts memory usage and generation time proportionally with negligible effect on the overdensity estimate. Default is 10.

• n_jobs (int) – Number of parallel worker processes for the (embarrassingly parallel, GLASS/healpy-bound) mock loop. 1 (default) runs serially; -1 uses all cores. Mock k always uses seed + k, so results are reproducible regardless of n_jobs.

Returns:

result –

• "delta_chi2" — shape (n_sys,), Δχ² for the data.

• "p_values" — shape (n_sys,), mock-based p-values in (0, 1]. Small values indicate significant systematic contamination.

• "delta_chi2_mocks" — shape (n_mocks, n_sys), Δχ² distribution from systematic-free mocks.

Return type:

dict with keys

Notes

Mocks are generated with generate_glass_fullsky_mock(), pixelised with pixelize_catalog(), and restricted to the survey footprint via good_pixels. The same template maps delta_t are used for both the data and the mocks.

p-value formula (smoothed to avoid zero):

\[p_i = \frac{\#\{k : \Delta\chi^2_{\rm mock,k,i} \ge \Delta\chi^2_{\rm data,i}\} + 1}{n_{\rm mocks} + 1}\]

Examples

>>> import numpy as np
>>> from sys_mapping import isd_template_significance
>>> rng = np.random.default_rng(0)
>>> n_pix, n_sys = 3072, 3   # nside=16 footprint size
>>> delta_t = rng.standard_normal((n_sys, n_pix))
>>> delta_g = 0.4 * delta_t[1] + rng.standard_normal(n_pix) * 0.1
>>> good = np.ones(12 * 16 ** 2, dtype=bool)
>>> z_edges = np.array([0.0, 0.5])
>>> nz = np.array([1000.0])
>>> result = isd_template_significance(
...     delta_g, delta_t, good, nside=16, n_total=3000,
...     z_edges=z_edges, nz=nz, n_mocks=5, seed=0)
>>> set(result.keys()) == {"delta_chi2", "p_values", "delta_chi2_mocks"}
True
>>> result["delta_chi2"].shape
(3,)
>>> result["p_values"].shape
(3,)
>>> result["delta_chi2_mocks"].shape
(5, 3)

References

Rodríguez-Monroy et al. 2025, arXiv:2509.07943, Sec. IV.A.1. Tessore et al. 2023, OJAp, 6, 11.