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)[source]

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')[source]

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

Three SNR definitions adapted from Tanidis et al. 2026 (Sec. 3) for galaxy clustering:

"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.

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".

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.

sys_mapping.diagnostics.footprint_mask_diagnostics(delta_g_obs, delta_t, mask_fractions, good_pixels)[source]

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.