sys_mapping.correction

Post-inference corrections: noise debiasing, template PCA rotation, parameter back-transformation, and two-point function correction.

The typical post-MCMC workflow is:

  1. debias_params() — subtract noise bias from squared MAP estimates.

  2. correct_two_point_function() — apply the debiased parameters to the observed \(w(\theta)\).

  3. correct_power_spectrum_harmonic() — apply the harmonic-space correction (Elsner et al. 2016).

Key papers: Berlfein et al. 2024 (Eq. 15–16, 21, Appendix A); Elsner et al. 2016 — see also Methods Reference.

Parameter debiasing, two-point function correction, and harmonic-space correction.

Implements:

  • Noise debiasing (Eq. 21): â²_i_debiased = â²_i - Var[â_i]

  • Template rotation (Appendix A): PCA diagonalizes template covariance

  • Two-point correction (Eq. 15-16)

  • Harmonic-space power spectrum correction (Elsner+2016)

sys_mapping.correction.debias_params(a_hat, b_hat, var_a, var_b)[source]

Debias squared contamination parameters (Eq. 21).

Corrects for noise inflation: E[â²] = a² + Var[â], so a²_debiased = â² − Var[â]. Negative values (when noise dominates) are clipped to zero (physically, squared amplitudes ≥ 0).

Parameters:

  • a_hat ((n_sys,) MLE / posterior median additive parameters)

  • b_hat ((n_sys,) MLE / posterior median multiplicative parameters)

  • var_a ((n_sys,) posterior variance of a_hat)

  • var_b ((n_sys,) posterior variance of b_hat)

Returns:

  • a_sq_debiased ((n_sys,) clipped to [0, ∞))

  • b_sq_debiased ((n_sys,) clipped to [0, ∞))

  • Performance

  • ———–

  • Measured on CPU (n_sys=5) (*~4 μs/call. Pure NumPy element-wise ops.*)

  • Precision

  • ———

  • When the exact posterior variance is known, the debiased estimate recovers

  • the true a² to atol < 1e-10 (limited by float64 round-off).

Return type:

tuple[ndarray, ndarray]

Examples

>>> import numpy as np
>>> from sys_mapping import debias_params
>>> a_hat = np.array([ 0.10, -0.05,  0.03])
>>> b_hat = np.array([ 0.08,  0.02, -0.04])
>>> var_a = np.array([0.002, 0.001, 0.0005])
>>> var_b = np.array([0.001, 0.0005, 0.002])
>>> a_sq, b_sq = debias_params(a_hat, b_hat, var_a, var_b)
>>> a_sq  # â² - Var[â], clipped at 0
array([0.008 , 0.0015, 0.0004])
>>> b_sq
array([0.005, 0.    , 0.    ])
sys_mapping.correction.rotate_templates(delta_t)[source]

Rotate templates to their PCA eigenbasis (Appendix A).

When templates are correlated, the likelihood parameters are degenerate. PCA rotation produces orthogonal templates so each parameter is identifiable. Eigenvectors are sorted by descending eigenvalue.

Parameters:

delta_t ((n_sys, n_pix) template values at unmasked pixels)

Returns:

  • delta_t_rot ((n_sys, n_pix) rotated (orthogonal) templates)

  • rotation_matrix ((n_sys, n_sys) V^T — rows are eigenvectors)

  • eigenvalues ((n_sys,) template variances in the rotated basis (descending))

  • Performance

  • ———–

  • Measured on CPU (n_sys=5, n_pix=10_000) (*~177 μs/call.*)

  • Uses numpy.linalg.eigh; scales as O(n_sys³ + n_sys × n_pix).

  • Precision

  • ———

  • Rotated template covariance is diagonal to atol < 1e-10.

  • Rotation matrix is orthogonal (R @ R.T = I) to atol < 1e-12.

  • Total variance (sum(eigenvalues)) equals trace(cov) to rtol < 1e-10.

Return type:

tuple[ndarray, ndarray, ndarray]

Examples

>>> import numpy as np
>>> from sys_mapping import rotate_templates
>>> rng = np.random.default_rng(0)
>>> delta_t = rng.standard_normal((3, 5000))
>>> delta_t_rot, R, eigvals = rotate_templates(delta_t)
>>> delta_t_rot.shape
(3, 5000)
>>> eigvals.shape   # sorted descending
(3,)
>>> # Verify orthogonality of rotated templates
>>> cov_rot = delta_t_rot @ delta_t_rot.T / 5000
>>> np.max(np.abs(cov_rot - np.diag(np.diag(cov_rot)))) < 1e-10
True
sys_mapping.correction.transform_params_from_rotated(a_rot, b_rot, rotation_matrix)[source]

Transform inferred parameters from rotated back to original template basis.

Applies a_orig = V @ a_rot where V = rotation_matrix.T (since rotation_matrix = V^T).

Parameters:

  • a_rot ((n_sys,) additive parameters in the rotated PCA basis)

  • b_rot ((n_sys,) multiplicative parameters in the rotated PCA basis)

  • rotation_matrix ((n_sys, n_sys) returned by rotate_templates())

Returns:

  • a_orig ((n_sys,) additive parameters in the original template basis)

  • b_orig ((n_sys,) multiplicative parameters in the original template basis)

  • Performance

  • ———–

  • Measured on CPU (n_sys=5) (*~4 μs/call. Two small matrix–vector products.*)

  • Precision

  • ———

  • Roundtrip (rotate → infer → back-transform) recovers original parameters

  • to atol < 1e-12 (float64 matrix–vector product).

Return type:

tuple[ndarray, ndarray]

Examples

>>> import numpy as np
>>> from sys_mapping import rotate_templates, transform_params_from_rotated
>>> rng = np.random.default_rng(0)
>>> delta_t = rng.standard_normal((3, 5000))
>>> _, R, _ = rotate_templates(delta_t)
>>> a = np.array([0.05, -0.03, 0.02])
>>> b = np.array([0.04,  0.01, -0.02])
>>> a_rot, b_rot = R @ a, R @ b
>>> a_back, b_back = transform_params_from_rotated(a_rot, b_rot, R)
>>> np.allclose(a_back, a, atol=1e-12)
True
sys_mapping.correction.correct_two_point_function(w_obs, a_hat, b_hat, var_a, var_b, template_correlations)[source]

Full pipeline: debias parameters, then correct two-point function.

Convenience wrapper that chains debias_params() → compute_two_point_correction().

Parameters:

  • w_obs ((n_bins,) observed angular correlation function)

  • a_hat ((n_sys,) posterior median additive parameters)

  • b_hat ((n_sys,) posterior median multiplicative parameters)

  • var_a ((n_sys,) posterior variance of a_hat (from get_param_variance_from_chain()))

  • var_b ((n_sys,) posterior variance of b_hat)

  • template_correlations ((n_sys, n_bins) ⟨δ_{t,i}(θ) δ_{t,i}(θ)⟩ per template per bin)

Returns:

  • (n_bins,) corrected angular correlation function

  • Performance

  • ———–

  • Measured on CPU (n_sys=5, n_bins=15) (*~1084 μs/call (JAX dispatch + debias).*)

  • Precision

  • ———

  • Propagates float64 NumPy arithmetic; limited by the precision of the

  • input estimates a_hat, b_hat, var_a, var_b.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from sys_mapping import correct_two_point_function
>>> n_sys, n_bins = 3, 15
>>> rng = np.random.default_rng(0)
>>> w_obs  = rng.standard_normal(n_bins) * 0.01
>>> a_hat  = np.array([0.05, -0.03, 0.02])
>>> b_hat  = np.array([0.04,  0.01, -0.02])
>>> var_a  = np.full(n_sys, 5e-4)
>>> var_b  = np.full(n_sys, 3e-4)
>>> tcorr  = np.abs(rng.standard_normal((n_sys, n_bins))) * 1e-3
>>> w_corr = correct_two_point_function(w_obs, a_hat, b_hat, var_a, var_b, tcorr)
>>> w_corr.shape
(15,)
sys_mapping.correction.correct_power_spectrum_harmonic(pseudo_cl, ell, n_templates, alpha, template_cls)[source]

Full harmonic-space correction pipeline (Elsner+2016).

Subtracts template pseudo-Cℓ contributions and applies the closed-form harmonic bias correction from template subtraction.

The three-step pipeline:

  1. Template subtraction: \(\tilde C_\ell^{\rm TS} = \hat C_\ell - \sum_i \hat\alpha_i\,C_\ell^{t_i}\).

  2. Harmonic bias correction: subtract \(b_\ell = -n/(2\ell+1)\).

Parameters:

  • pseudo_cl (ndarray) – Measured pseudo-Cℓ of the galaxy field (shape (n_ell,)).

  • ell (ndarray) – Multipole array (shape (n_ell,)).

  • n_templates (int) – Number of subtracted templates (must equal len(alpha)).

  • alpha (ndarray) – Contamination amplitude estimates (shape (n_templates,)).

  • template_cls (ndarray) – Per-template pseudo-Cℓ (shape (n_templates, n_ell)). template_cls[i] is the pseudo-Cℓ of template map i.

Returns:

Bias-corrected power spectrum estimate (shape (n_ell,)).

Return type:

cl_corrected

Notes

The template pseudo-Cℓ can be pre-computed from the masked template maps via measure_pseudo_cl().

Examples

>>> import numpy as np
>>> from sys_mapping import correct_power_spectrum_harmonic
>>> n_ell, n_sys = 30, 2
>>> ell = np.arange(n_ell)
>>> cl_obs = np.ones(n_ell) * 1e-4
>>> alpha = np.array([0.1, -0.05])
>>> t_cls = np.ones((n_sys, n_ell)) * 5e-5
>>> cl_corr = correct_power_spectrum_harmonic(cl_obs, ell, n_sys, alpha, t_cls)
>>> cl_corr.shape
(30,)

References

Elsner, Leistedt & Peiris 2016, MNRAS 456, 2095. Elsner, Leistedt & Peiris 2017, MNRAS 465, 1847.