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
.
$ \tilde{I} = \int g(\nu) I(\nu)d\nu$
If we consider emission of the CMB in $K_{CMB}$ units, $I_{CMB}$ is not a function of $\nu$:
$ \tilde{I}_{CMB}[K_{CMB}] = \int g(\nu) I_{CMB}[K_{CMB}] d\nu = I_{CMB}[K_{CMB}] \int g(\nu) d\nu = I_{CMB}[K_{CMB}]$
$ \tilde{I}_{CMB}[K_{RJ}] = I_{CMB}[K_{CMB}] \int g(\nu) C_{K_{CMB}}^{K_{RJ}}(\nu) d\nu $
However we assume that the bandpass of the instrument $g(\nu)$ is given in power units, i.e. $Jy/sr$, the python normalize_weights
functions does:
$ g [K_{RJ}/K_{RJ}] = \dfrac{C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}]}{\int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}] d(\nu)} $
CMB bandpass integrated converted to $K_{RJ}$.
$ \tilde{I}_{CMB}[K_{RJ}] = I_{CMB}[K_{CMB}] \int g(\nu)[K_{RJ}/K_{RJ}] C_{K_{CMB}}^{K_{RJ}}(\nu) d\nu = I_{CMB}[K_{CMB}] \int \dfrac {C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}]}{\int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}] d(\nu)} C_{K_{CMB}}^{K_{RJ}}(\nu) d\nu $
You can think the last equation as getting a value in $K_{CMB}$, turn it into $K_{RJ}$ and then to $Jy~sr^{-1}$, and do the relative weighting then.
However we assume that the bandpass of the instrument $g(\nu)$ is given in power units, i.e. $Jy/sr$, the python normalize_weights
functions does:
$ g [K_{RJ}/K_{RJ}] = \dfrac{C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}]}{\int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}] d(\nu)} $
If the output is requested in $K_{CMB}$, we basically have to undo that steps, so:
$ \tilde{I}_{CMB}[K_{RJ}] = I_{CMB}[K_{CMB}] \int \dfrac{C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}]}{\int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}] d(\nu)} C_{K_{CMB}}^{K_{RJ}}(\nu) d\nu $
fixed_bandpass_unit_conversion
¶$ \tilde{I}_{CMB}[K_{CMB}] = \tilde{I}_{CMB}[K_{RJ}] \dfrac{ \int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu} { \int C_{K_{CMB}}^{Jy~sr^{-1}}(\nu) g(\nu) d\nu} $
broken_bandpass_unit_conversion
¶$ \tilde{I}_{CMB}[K_{CMB}] = \tilde{I}_{CMB}[K_{RJ}] \int C_{K_{RJ}}^{K_{CMB}}(\nu) g_{RJ}(\nu) d(\nu) = \tilde{I}_{CMB}[K_{RJ}] \int C_{K_{RJ}}^{K_{CMB}}(\nu) \dfrac{C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}]}{\int C_{K_{RJ}}^{Jy~sr^{-1}}(\nu) g [Jy~sr^{-1} / Jy~sr^{-1}] d(\nu)} d(\nu)$
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)
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%
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%
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%