Investigate broken implementation of bandpass unit conversion in PySM 3

import pysm3
import numpy as np
import healpy as hp
from astropy.tests.helper import assert_quantity_allclose
import pysm3.units as u

Investigation of the broken implementation of the bandpass unit conversion function in PySM 3 and check of its impact on existing Simons Observatory and CMB-S4 simulated datasets.

See the related issue and the pull request with the fix. This affects all PySM 3 versions <3.2.2, it will be fixed in 3.3.0.

$ = g() I()d$

If we consider emission of the CMB in $K_{C M B}$ units, $I_{C M B}$ is not a function of $ν$ :

$ {CMB}[K_{CMB}] = g() I{CMB}[K_{CMB}] d= I_{CMB}[K_{CMB}] g() d= I_{CMB}[K_{CMB}]$

$ {CMB}[K_{RJ}] = I{CMB}[K_{CMB}] g() C_{K_{CMB}}^{K_{RJ}}() d$

However we assume that the bandpass of the instrument $g (ν)$ is given in power units, i.e. $J y / s r$ , the python normalize_weights functions does:

$ g [K_{RJ}/K_{RJ}] = $

CMB bandpass integrated converted to $K_{R J}$ .

$ {CMB}[K_{RJ}] = I{CMB}[K_{CMB}] g()[K_{RJ}/K_{RJ}] C_{K_{CMB}}^{K_{RJ}}() d= I_{CMB}[K_{CMB}] C_{K_{CMB}}^{K_{RJ}}() d$

You can think the last equation as getting a value in $K_{C M B}$ , turn it into $K_{R J}$ and then to $J y s r^{- 1}$ , and do the relative weighting then.

However we assume that the bandpass of the instrument $g (ν)$ is given in power units, i.e. $J y / s r$ , the python normalize_weights functions does:

$ g [K_{RJ}/K_{RJ}] = $

From $K_{C M B}$ to $K_{R J}$ :

If the output is requested in $K_{C M B}$ , we basically have to undo that steps, so:

$ {CMB}[K_{RJ}] = I{CMB}[K_{CMB}] C_{K_{CMB}}^{K_{RJ}}() d$

`fixed_bandpass_unit_conversion`

$ {CMB}[K_{CMB}] = {CMB}[K_{RJ}] { C_{K_{CMB}}^{Jy~sr{-1}}() g() d} $

`broken_bandpass_unit_conversion`

$ {CMB}[K_{CMB}] = {CMB}[K_{RJ}] C_{K_{RJ}}^{K_{CMB}}() g_{RJ}() d() = {CMB}[K_{RJ}] C{K_{RJ}}^{K_{CMB}}() d()$

from pysm3.utils import normalize_weights

def broken_bandpass_unit_conversion(freqs, weights, output_unit):
    """Unit conversion from uK_RJ to output unit given a bandpass

    Parameters
    ----------
    freqs : astropy.units.Quantity
        Frequency array in a unit compatible with GHz
    """
    freqs = check_freq_input(freqs)
    factors = (np.ones(len(freqs), dtype=np.float) * u.uK_RJ).to_value(
        output_unit, equivalencies=u.cmb_equivalencies(freqs * u.GHz)
    )
    if len(freqs) > 1:
        w = normalize_weights(freqs, weights)
        factor = np.trapz(factors * w, freqs)
    else:
        factor = factors[0]
    return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)

    nside = 32
    freqs = np.array([250, 300, 350]) * u.GHz
    weights = np.ones(len(freqs))
    sky = pysm3.Sky(nside=nside, preset_strings=["c2"])
    CMB_rj_int = sky.get_emission(freqs, weights)
    CMB_thermo_int = CMB_rj_int*fixed_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    expected_map = pysm3.read_map(
        "pysm_2/lensed_cmb.fits", field=(0, 1), nside=nside, unit=u.uK_CMB
    )
    for pol in [0, 1]:
        #assert_quantity_allclose(expected_map[pol], CMB_thermo_int[pol], rtol=1e-4)
        print("Error: {:.2f} %".format(100-np.median((CMB_thermo_int[pol]/expected_map[pol]).value)*100))

Error: -0.00 %
Error: -0.00 %

from pysm3.utils import check_freq_input

def fixed_bandpass_unit_conversion(freqs, weights, output_unit):
    """Unit conversion from uK_RJ to output unit given a bandpass
    Parameters
    ----------
    freqs : astropy.units.Quantity
        Frequency array in a unit compatible with GHz
    """
    freqs = check_freq_input(freqs)
    weights_to_rj = (weights * u.uK_RJ).to_value(
            (u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
        )
    weights_to_out = (weights * output_unit).to_value(
            (u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
        )
    if len(freqs) > 1:
        factor = np.trapz(weights_to_rj, freqs)/np.trapz(weights_to_out, freqs)
    else:
        factor = (1. * u.uK_RJ).to_value(
            output_unit, equivalencies=u.cmb_equivalencies(freqs * u.GHz)
        )
    return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)

Unit conversion error for 20% tophat bandpasses

perc_error = {}
for center_freq in np.array([10, 20, 40, 70, 100, 150, 200, 250, 270, 300, 350, 400])*u.GHz:
    freqs = np.linspace(center_freq*.8, center_freq*1.2, 10)
    weights = np.ones(10)
    rj_to_cmb = broken_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
    print("{: <3.0f}, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(center_freq, rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))

10  GHz, broken: 1.003 uK_CMB / uK_RJ, fixed: 1.003 uK_CMB / uK_RJ, error 0.00%
20  GHz, broken: 1.011 uK_CMB / uK_RJ, fixed: 1.011 uK_CMB / uK_RJ, error 0.00%
40  GHz, broken: 1.045 uK_CMB / uK_RJ, fixed: 1.045 uK_CMB / uK_RJ, error 0.01%
70  GHz, broken: 1.143 uK_CMB / uK_RJ, fixed: 1.142 uK_CMB / uK_RJ, error 0.08%
100 GHz, broken: 1.310 uK_CMB / uK_RJ, fixed: 1.306 uK_CMB / uK_RJ, error 0.32%
150 GHz, broken: 1.808 uK_CMB / uK_RJ, fixed: 1.782 uK_CMB / uK_RJ, error 1.47%
200 GHz, broken: 2.770 uK_CMB / uK_RJ, fixed: 2.661 uK_CMB / uK_RJ, error 4.06%
250 GHz, broken: 4.627 uK_CMB / uK_RJ, fixed: 4.260 uK_CMB / uK_RJ, error 8.62%
270 GHz, broken: 5.800 uK_CMB / uK_RJ, fixed: 5.221 uK_CMB / uK_RJ, error 11.10%
300 GHz, broken: 8.297 uK_CMB / uK_RJ, fixed: 7.178 uK_CMB / uK_RJ, error 15.59%
350 GHz, broken: 15.749 uK_CMB / uK_RJ, fixed: 12.560 uK_CMB / uK_RJ, error 25.39%
400 GHz, broken: 31.306 uK_CMB / uK_RJ, fixed: 22.588 uK_CMB / uK_RJ, error 38.59%

Unit conversion error for Simons Observatory channels

import h5py

SO_chs = h5py.File("../mapsims/mapsims/data/simonsobs_instrument_parameters_2020.06.h5", mode='r')

SO_chs["LT0_UHF1"].keys()

<KeysViewHDF5 ['bandpass_frequency_GHz', 'bandpass_weight']>

SO_chs["LT0_UHF1"].attrs.keys()

<KeysViewHDF5 ['band', 'band_id', 'center_frequency_GHz', 'fwhm_arcmin', 'noise_band_index', 'telescope', 'tube', 'tube_id']>

perc_error = {}
for ch in SO_chs.values():
    freqs = ch["bandpass_frequency_GHz"] * u.GHz
    weights = ch["bandpass_weight"]
    rj_to_cmb = broken_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
    print("{} {:4s}, {:6.2f} GHz, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(ch.attrs["telescope"], ch.attrs["band"], ch.attrs["center_frequency_GHz"], rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))

LA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
LA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
LA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
LA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
LA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
LA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
LA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
LA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
LA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
LA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
LA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
LA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
LA LF1 ,  25.70 GHz, broken: 1.018 uK_CMB / uK_RJ, fixed: 1.018 uK_CMB / uK_RJ, error 0.00%
LA LF2 ,  38.90 GHz, broken: 1.043 uK_CMB / uK_RJ, fixed: 1.043 uK_CMB / uK_RJ, error 0.01%
SA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
SA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
SA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
SA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
SA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
SA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
SA LF1 ,  25.70 GHz, broken: 1.018 uK_CMB / uK_RJ, fixed: 1.018 uK_CMB / uK_RJ, error 0.00%
SA LF2 ,  38.90 GHz, broken: 1.043 uK_CMB / uK_RJ, fixed: 1.043 uK_CMB / uK_RJ, error 0.01%

Unit conversion error for CMB-S4 channels

S4_chs = h5py.File("../s4mapbasedsims/202006_foregrounds_extragalactic_cmb_tophat/cmbs4_tophat.h5", mode='r')

perc_error = {}
for ch in S4_chs.values():
    freqs = ch["bandpass_frequency_GHz"] * u.GHz
    weights = ch["bandpass_weight"]
    rj_to_cmb = broken_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
        freqs, weights, u.uK_CMB
    )
    perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
    print("{} {:5s}, {:6.2f} GHz, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(
        ch.attrs["telescope"], ch.attrs["band"], ch.attrs["center_frequency_GHz"], rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))

LAT HFL1 , 225.00 GHz, broken: 3.364 uK_CMB / uK_RJ, fixed: 3.275 uK_CMB / uK_RJ, error 2.71%
LAT HFL2 , 278.00 GHz, broken: 5.632 uK_CMB / uK_RJ, fixed: 5.523 uK_CMB / uK_RJ, error 1.98%
SAT HFS1 , 220.00 GHz, broken: 3.163 uK_CMB / uK_RJ, fixed: 3.109 uK_CMB / uK_RJ, error 1.71%
SAT HFS2 , 270.00 GHz, broken: 5.269 uK_CMB / uK_RJ, fixed: 5.099 uK_CMB / uK_RJ, error 3.34%
LAT LFL1 ,  27.00 GHz, broken: 1.019 uK_CMB / uK_RJ, fixed: 1.019 uK_CMB / uK_RJ, error 0.00%
LAT LFL2 ,  39.00 GHz, broken: 1.044 uK_CMB / uK_RJ, fixed: 1.044 uK_CMB / uK_RJ, error 0.01%
SAT LFS1 ,  30.00 GHz, broken: 1.024 uK_CMB / uK_RJ, fixed: 1.024 uK_CMB / uK_RJ, error 0.00%
SAT LFS2 ,  40.00 GHz, broken: 1.044 uK_CMB / uK_RJ, fixed: 1.044 uK_CMB / uK_RJ, error 0.01%
SAT MFHS1,  95.00 GHz, broken: 1.263 uK_CMB / uK_RJ, fixed: 1.262 uK_CMB / uK_RJ, error 0.10%
SAT MFHS2, 155.10 GHz, broken: 1.822 uK_CMB / uK_RJ, fixed: 1.813 uK_CMB / uK_RJ, error 0.51%
LAT MFL1 ,  93.00 GHz, broken: 1.262 uK_CMB / uK_RJ, fixed: 1.259 uK_CMB / uK_RJ, error 0.22%
LAT MFL2 , 145.00 GHz, broken: 1.707 uK_CMB / uK_RJ, fixed: 1.697 uK_CMB / uK_RJ, error 0.62%
SAT MFLS1,  85.00 GHz, broken: 1.207 uK_CMB / uK_RJ, fixed: 1.206 uK_CMB / uK_RJ, error 0.06%
SAT MFLS2, 145.10 GHz, broken: 1.696 uK_CMB / uK_RJ, fixed: 1.690 uK_CMB / uK_RJ, error 0.40%
LAT ULFL1,  20.00 GHz, broken: 1.011 uK_CMB / uK_RJ, fixed: 1.011 uK_CMB / uK_RJ, error 0.00%

From KCMB to KRJ: