From 74c815d80679d4820208b13b69f0fc80e83e782e Mon Sep 17 00:00:00 2001 From: brandonshensley Date: Mon, 4 Oct 2021 20:47:44 -0400 Subject: [PATCH 01/44] Fixed year on Bennett et al. citation --- docs/models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/models.rst b/docs/models.rst index 42292339..ca5a5e34 100644 --- a/docs/models.rst +++ b/docs/models.rst @@ -53,7 +53,7 @@ Dust Synchrotron =========== -- **s1**: A power law scaling is used for the synchrotron emission, with a spatially varying spectral index. The emission templates are the Haslam 408 MHz, 57' resolution data reprocessed by Remazeilles et al 2015 MNRAS 451, 4311, and the WMAP 9-year 23 GHz Q/U maps (Bennett, C.L., et.al., 2014, ApJS, 208, 20B). The polarization maps have been smoothed with a Gaussian kernel of FWHM 5 degrees and had small scales added. The intensity template has had small scales added straight to the template. The details of the small scale procedure is outlined in the accompanying paper. The spectral index map was derived using a combination of the Haslam 408 MHz data and WMAP 23 GHz 7-year data (Miville-Deschenes, M.-A. et al., 2008, A&A, 490, 1093). The same scaling is used for intensity and polarization. This is the same prescription as used in the Planck Sky Model's v1.7.8 'power law' option (Delabrouille et al. A&A 553, A96, 2013), but with the Haslam map updated to the Remazeilles version. A 'curved power law' model is also supported with a single isotropic curvature index. The amplitude of this curvature is taken from Kogut, A. 2012, ApJ, 753, 110. +- **s1**: A power law scaling is used for the synchrotron emission, with a spatially varying spectral index. The emission templates are the Haslam 408 MHz, 57' resolution data reprocessed by Remazeilles et al 2015 MNRAS 451, 4311, and the WMAP 9-year 23 GHz Q/U maps (Bennett, C.L., et.al., 2013, ApJS, 208, 20B). The polarization maps have been smoothed with a Gaussian kernel of FWHM 5 degrees and had small scales added. The intensity template has had small scales added straight to the template. The details of the small scale procedure is outlined in the accompanying paper. The spectral index map was derived using a combination of the Haslam 408 MHz data and WMAP 23 GHz 7-year data (Miville-Deschenes, M.-A. et al., 2008, A&A, 490, 1093). The same scaling is used for intensity and polarization. This is the same prescription as used in the Planck Sky Model's v1.7.8 'power law' option (Delabrouille et al. A&A 553, A96, 2013), but with the Haslam map updated to the Remazeilles version. A 'curved power law' model is also supported with a single isotropic curvature index. The amplitude of this curvature is taken from Kogut, A. 2012, ApJ, 753, 110. - **s2**: synchrotron index steepens off the Galactic plane, from -3.0 in the plane to -3.3 off the plane. Consistent with WMAP. From ec529ccc59ff0996f011212ec401e0d18305347f Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 17:29:50 -0700 Subject: [PATCH 02/44] Remove doctest-rst option from pytest We do not have any code block in the documentation See https://github.com/openjournals/joss-reviews/issues/3783#issuecomment-932672403 --- setup.cfg | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index d147c94f..aba0981f 100644 --- a/setup.cfg +++ b/setup.cfg @@ -42,7 +42,6 @@ testpaths = "pysm3" "docs" astropy_header = true doctest_plus = enabled text_file_format = rst -addopts = --doctest-rst [coverage:run] omit = From f578230187c173aae9f67be8afa96d6540bd6908 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 18:34:20 -0700 Subject: [PATCH 03/44] Rewrite test_units.py Use pytest instead of unittest, simplify using tmp_path --- pysm3/tests/test_units.py | 151 ++++++++++++++++---------------------- 1 file changed, 63 insertions(+), 88 deletions(-) diff --git a/pysm3/tests/test_units.py b/pysm3/tests/test_units.py index 4e80e644..3b8212fe 100644 --- a/pysm3/tests/test_units.py +++ b/pysm3/tests/test_units.py @@ -3,105 +3,80 @@ import numpy as np import healpy as hp import pysm3.units as units +import pytest from pysm3.models.template import read_map -class TestUnits(unittest.TestCase): - def setUp(self): - self.T_CMB = 100.0 * units.K_CMB - self.T_RJ = 100.0 * units.K_RJ - self.freqs = 100.0 * units.GHz - self.freqs_arr = np.array([10.0, 100.0, 1000.0]) * units.GHz - self.temp_dir = Path(__file__).absolute().parent / "temp_dir" - self.temp_dir.mkdir(parents=True, exist_ok=True) - self.temp_fits_file_RJ = self.temp_dir / "test_RJ.fits" - self.temp_fits_file_CMB = self.temp_dir / "test_CMB.fits" - self.temp_fits_file_dimless = self.temp_dir / "test_dimless.fits" - self.temp_fits_file_no_unit_hdr = self.temp_dir / "test_no_hdr.fits" - nside = 256 - npix = hp.nside2npix(nside) - self.test_map_RJ = np.random.randn(npix) * units.K_RJ - self.test_map_CMB = np.random.randn(npix) * units.K_CMB - self.test_map_dimless = units.Quantity(np.random.randn(npix), "") +@pytest.fixture +def setUp(): + T_CMB = 100.0 * units.K_CMB + T_RJ = 100.0 * units.K_RJ + freqs = 100.0 * units.GHz + nside = 256 + npix = hp.nside2npix(nside) + test_map_RJ = np.random.randn(npix) * units.K_RJ + test_map_CMB = np.random.randn(npix) * units.K_CMB + test_map_dimless = units.Quantity(np.random.randn(npix), "") + return T_CMB, T_RJ, freqs, test_map_RJ, test_map_CMB, test_map_dimless - def tearDown(self): - try: - self.temp_fits_file_RJ.unlink() - except FileNotFoundError: - pass - try: - self.temp_fits_file_CMB.unlink() - except FileNotFoundError: - pass - try: - self.temp_fits_file_dimless.unlink() - except FileNotFoundError: - pass - try: - self.temp_fits_file_no_unit_hdr.unlink() - except FileNotFoundError: - pass - self.temp_dir.rmdir() - def test_conversion(self): - """ Here we test that the numerical value of the conversion is correct. - The mathematical form is +def test_conversion(setUp): + """ Here we test that the numerical value of the conversion is correct. + The mathematical form is - ..math:: - I_\nu = \frac{2 \nu^2 k T_{\rm RJ}}{c^2} \\ - I_\nu = T_{\rm CMB} B^\prime_\nu(T_{\rm CMB, 0}) + ..math:: + I_\nu = \frac{2 \nu^2 k T_{\rm RJ}}{c^2} \\ + I_\nu = T_{\rm CMB} B^\prime_\nu(T_{\rm CMB, 0}) - so, eliminating the flux in this equation: + so, eliminating the flux in this equation: - ..math:: - T_{\rm RJ} / T_{\rm CMB} = \frac{c^2}{2 \nu^2 k_B}B^\prime_\nu(T_{\rm CMB, 0}) + ..math:: + T_{\rm RJ} / T_{\rm CMB} = \frac{c^2}{2 \nu^2 k_B}B^\prime_\nu(T_{\rm CMB, 0}) - Here we calculate the RHS of this equation and compare it to the - ratio of T_RJ and the result of its transformation to T_CMB. - """ - equiv = {"equivalencies": units.cmb_equivalencies(self.freqs)} - rj_from_cmb = self.T_CMB.to(units.K_RJ, **equiv) - cmb_from_rj = self.T_RJ.to(units.K_CMB, **equiv) + Here we calculate the RHS of this equation and compare it to the + ratio of T_RJ and the result of its transformation to T_CMB. + """ + T_CMB, T_RJ, freqs, test_map_RJ, test_map_CMB, test_map_dimless = setUp + equiv = {"equivalencies": units.cmb_equivalencies(freqs)} + rj_from_cmb = T_CMB.to(units.K_RJ, **equiv) + cmb_from_rj = T_RJ.to(units.K_CMB, **equiv) - # check that the reverse transformation gives overall transformation of unity. - reverse1 = rj_from_cmb.to(units.K_CMB, **equiv) - reverse2 = cmb_from_rj.to(units.K_RJ, **equiv) + # check that the reverse transformation gives overall transformation of unity. + reverse1 = rj_from_cmb.to(units.K_CMB, **equiv) + reverse2 = cmb_from_rj.to(units.K_RJ, **equiv) - np.testing.assert_almost_equal(1.0, self.T_CMB / reverse1, decimal=6) - np.testing.assert_almost_equal(1.0, self.T_RJ / reverse2, decimal=6) + np.testing.assert_almost_equal(1.0, T_CMB / reverse1, decimal=6) + np.testing.assert_almost_equal(1.0, T_RJ / reverse2, decimal=6) - def test_fits_unit_funcitonality(self): - """ Test that the units can be written to the fits header. Check that - they can be read in again and assigned to the data in that fits file - correctly. - """ - hp.write_map( - str(self.temp_fits_file_RJ), - self.test_map_RJ.value, - column_units=self.test_map_RJ.unit.to_string("generic"), - ) - hp.write_map( - str(self.temp_fits_file_CMB), - self.test_map_CMB.value, - column_units=self.test_map_CMB.unit.to_string("generic"), - ) - hp.write_map( - str(self.temp_fits_file_dimless), - self.test_map_dimless.value, - column_units=self.test_map_dimless.unit.to_string("generic"), - ) - hp.write_map(str(self.temp_fits_file_no_unit_hdr), self.test_map_dimless.value) - cmb_in = read_map(str(self.temp_fits_file_CMB), 256) - rj_in = read_map(str(self.temp_fits_file_RJ), 256) - dimless_in = read_map(str(self.temp_fits_file_dimless), 256) - no_unit_hdr = read_map(str(self.temp_fits_file_no_unit_hdr), 256) - self.assertTrue(cmb_in.unit == units.K_CMB) - self.assertTrue(rj_in.unit == units.K_RJ) - self.assertTrue(dimless_in.unit == units.dimensionless_unscaled) - self.assertTrue(no_unit_hdr.unit == units.dimensionless_unscaled) - return +def test_fits_unit_functionality(setUp, tmp_path): + """Test that the units can be written to the fits header. Check that + they can be read in again and assigned to the data in that fits file + correctly. + """ + T_CMB, T_RJ, freqs, test_map_RJ, test_map_CMB, test_map_dimless = setUp + hp.write_map( + tmp_path / "temp_fits_file_RJ.fits", + test_map_RJ.value, + column_units=test_map_RJ.unit.to_string("generic"), + ) + hp.write_map( + tmp_path / "temp_fits_file_CMB.fits", + test_map_CMB.value, + column_units=test_map_CMB.unit.to_string("generic"), + ) + hp.write_map( + tmp_path / "temp_fits_file_dimless.fits", + test_map_dimless.value, + column_units=test_map_dimless.unit.to_string("generic"), + ) + hp.write_map(tmp_path / "temp_fits_file_no_unit_hdr.fits", test_map_dimless.value) - -if __name__ == "__main__": - unittest.main() + cmb_in = read_map(tmp_path / "temp_fits_file_CMB.fits", 256) + rj_in = read_map(tmp_path / "temp_fits_file_RJ.fits", 256) + dimless_in = read_map(tmp_path / "temp_fits_file_dimless.fits", 256) + no_unit_hdr = read_map(tmp_path / "temp_fits_file_no_unit_hdr.fits", 256) + assert cmb_in.unit == units.K_CMB + assert rj_in.unit == units.K_RJ + assert dimless_in.unit == units.dimensionless_unscaled + assert no_unit_hdr.unit == units.dimensionless_unscaled From 0139f1c545a558b39377c96e88ab404688cf50e2 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:50:58 -0700 Subject: [PATCH 04/44] fix: remove unused imports --- pysm3/tests/test_units.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pysm3/tests/test_units.py b/pysm3/tests/test_units.py index 3b8212fe..98a79cd9 100644 --- a/pysm3/tests/test_units.py +++ b/pysm3/tests/test_units.py @@ -1,5 +1,3 @@ -from pathlib import Path -import unittest import numpy as np import healpy as hp import pysm3.units as units From f4a8aa09b772259c5e74ba16c10ddc90cfb6711d Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 17:38:22 -0700 Subject: [PATCH 05/44] Use logging instead of warning --- pysm3/utils/data.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pysm3/utils/data.py b/pysm3/utils/data.py index 0665f0ab..b7f52b0f 100644 --- a/pysm3/utils/data.py +++ b/pysm3/utils/data.py @@ -1,5 +1,6 @@ import warnings import os +import logging from astropy.utils import data @@ -7,6 +8,7 @@ PREDEFINED_DATA_FOLDERS = ["/global/project/projectdirs/cmb/www/pysm-data/"] # NERSC +log = logging.getLogger("pysm3") class RemoteData: def __init__(self): @@ -34,11 +36,11 @@ def get(self, filename): for folder in self.data_folders: full_path = os.path.join(folder, filename) if os.path.exists(full_path): - warnings.warn(f"Access data from {full_path}") + log.info(f"Access data from {full_path}") return full_path with data.conf.set_temp("dataurl", self.data_url), data.conf.set_temp( "remote_timeout", 90 ): - warnings.warn(f"Retrieve data for {filename} (if not cached already)") + log.info(f"Retrieve data for {filename} (if not cached already)") full_path = data.get_pkg_data_filename(filename, show_progress=True) return full_path From 9cdafd823732c153fc5cf776e34f515caa44962c Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 17:40:39 -0700 Subject: [PATCH 06/44] Log instead of warning for unit --- pysm3/models/template.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pysm3/models/template.py b/pysm3/models/template.py index 94f0057e..3a40836b 100644 --- a/pysm3/models/template.py +++ b/pysm3/models/template.py @@ -7,7 +7,7 @@ Objects: Model """ -import warnings +import logging import numpy as np import healpy as hp from astropy.io import fits @@ -16,6 +16,7 @@ from .. import mpi import gc +log = logging.getLogger("pysm3") class Model: """ This is the template object for PySM objects. @@ -189,7 +190,7 @@ def extract_hdu_unit(path): except KeyError: # in the case that TUNIT1 does not exist, assume unitless quantity. unit = "" - warnings.warn("No physical unit associated with file " + str(path)) + log.warning("No physical unit associated with file %s", str(path)) return unit From ae0b83be40bfd46604f057920c58bc5ef5dc523f Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 18:38:56 -0700 Subject: [PATCH 07/44] fix codestyle --- pysm3/models/template.py | 9 +++++---- pysm3/utils/data.py | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/pysm3/models/template.py b/pysm3/models/template.py index 3a40836b..7cd47573 100644 --- a/pysm3/models/template.py +++ b/pysm3/models/template.py @@ -18,8 +18,9 @@ log = logging.getLogger("pysm3") + class Model: - """ This is the template object for PySM objects. + """This is the template object for PySM objects. If a MPI communicator is passed as input and `pixel_indices` is None, the class automatically distributes the maps across processes. @@ -68,7 +69,7 @@ def read_txt(self, path, **kwargs): @u.quantity_input def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ: - """ This function evaluates the component model at a either + """This function evaluates the component model at a either a single frequency, an array of frequencies, or over a bandpass. Parameters @@ -153,7 +154,7 @@ def apply_smoothing_and_coord_transform( def apply_normalization(freqs, weights): - """ Function to apply a normalization constraing to a set of weights. + """Function to apply a normalization constraing to a set of weights. This imposes the requirement that the integral of the weights over the array `freqs` must equal unity. @@ -174,7 +175,7 @@ def apply_normalization(freqs, weights): def extract_hdu_unit(path): - """ Function to extract unit from an hdu. + """Function to extract unit from an hdu. Parameters ---------- path: Path object diff --git a/pysm3/utils/data.py b/pysm3/utils/data.py index b7f52b0f..77779a99 100644 --- a/pysm3/utils/data.py +++ b/pysm3/utils/data.py @@ -1,4 +1,3 @@ -import warnings import os import logging @@ -10,6 +9,7 @@ log = logging.getLogger("pysm3") + class RemoteData: def __init__(self): """Access template from remote server From 6ca3741697705f378fda80bb85d56441d3f09b36 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:00:28 -0700 Subject: [PATCH 08/44] move also bandpass_unit_conversion to logging --- pysm3/utils/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pysm3/utils/__init__.py b/pysm3/utils/__init__.py index b102c1e9..dee80edb 100644 --- a/pysm3/utils/__init__.py +++ b/pysm3/utils/__init__.py @@ -3,14 +3,14 @@ # This sub-module is destined for common non-package specific utility # functions. -import warnings - import numpy as np from numba import njit from .. import units as u from .data import RemoteData # noqa: F401 +import logging +log = logging.getLogger("pysm3") def has_polarization(m): """Checks if a map or a group of map is polarized @@ -111,7 +111,7 @@ def bandpass_unit_conversion( weights /= np.trapz(weights, freqs) if weights.min() < cut: good = np.logical_not(weights < cut) - warnings.warn( + log.info( "Removing {}/{} points below {}".format(good.sum(), len(good), cut) ) weights = weights[good] From 799b886b58d28931028b295c3a29dd9efeb263be Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:04:00 -0700 Subject: [PATCH 09/44] docs: how to handle verbosity --- docs/index.rst | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/docs/index.rst b/docs/index.rst index 8d135cc3..f89766f5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -102,6 +102,33 @@ Execute the unit tests with:: pytest +Configure verbosity +------------------- + +PySM uses the `logging` module to configure its verbosity, +by default it will only print warnings and errors, to configure logging +you can access the "pysm3" logger with:: + + import logging + log = logging.getLogger("pysm3") + +configure the logging level:: + + log.setLevel(logging.DEBUG) + +redirect the logs to the console:: + + handler = logging.StreamHandler() + log.addHandler(handler) + +or customize their format:: + + log_format="%(name)s - %(levelname)s - %(message)s" + formatter = logging.Formatter(log_format) + handler.setFormatter(formatter) + +For more details see the `Python documentation `_. + Tutorials ========= From aec18c7b3d4346e69ebddb303d683a8902720d93 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:10:26 -0700 Subject: [PATCH 10/44] fix: escape \e that was causing DeprecationWarning --- pysm3/models/template.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysm3/models/template.py b/pysm3/models/template.py index 7cd47573..19767785 100644 --- a/pysm3/models/template.py +++ b/pysm3/models/template.py @@ -43,7 +43,7 @@ def __init__(self, nside, map_dist=None): nside: int Resolution parameter at which this model is to be calculated. smoothing_lmax : int - :math:`\ell_{max}` for the smoothing step, by default :math:`2*N_{side}` + :math:`\\ell_{max}` for the smoothing step, by default :math:`2*N_{side}` """ self.nside = nside assert nside is not None From f6eaaf020fb107cec60ecab2d4608a4a28dbd7e7 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:15:50 -0700 Subject: [PATCH 11/44] fix: extract_hdu_unit was leaving FITS file open --- pysm3/models/template.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pysm3/models/template.py b/pysm3/models/template.py index 19767785..006c511c 100644 --- a/pysm3/models/template.py +++ b/pysm3/models/template.py @@ -185,13 +185,13 @@ def extract_hdu_unit(path): string String specifying the unit of the fits data. """ - hdul = fits.open(path) - try: - unit = hdul[1].header["TUNIT1"] - except KeyError: - # in the case that TUNIT1 does not exist, assume unitless quantity. - unit = "" - log.warning("No physical unit associated with file %s", str(path)) + with fits.open(path) as hdul: + try: + unit = hdul[1].header["TUNIT1"] + except KeyError: + # in the case that TUNIT1 does not exist, assume unitless quantity. + unit = "" + log.warning("No physical unit associated with file %s", str(path)) return unit From 153a12cb2c63977014ec664af21cad64e0ebe9f6 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:24:58 -0700 Subject: [PATCH 12/44] escape back slash --- pysm3/tests/test_units.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pysm3/tests/test_units.py b/pysm3/tests/test_units.py index 98a79cd9..cf0b3f2f 100644 --- a/pysm3/tests/test_units.py +++ b/pysm3/tests/test_units.py @@ -23,13 +23,13 @@ def test_conversion(setUp): The mathematical form is ..math:: - I_\nu = \frac{2 \nu^2 k T_{\rm RJ}}{c^2} \\ - I_\nu = T_{\rm CMB} B^\prime_\nu(T_{\rm CMB, 0}) + I_\\nu = \\frac{2 \\nu^2 k T_{\\rm RJ}}{c^2} \\\\ + I_\\nu = T_{\\rm CMB} B^\\prime_\\nu(T_{\\rm CMB, 0}) so, eliminating the flux in this equation: ..math:: - T_{\rm RJ} / T_{\rm CMB} = \frac{c^2}{2 \nu^2 k_B}B^\prime_\nu(T_{\rm CMB, 0}) + T_{\\rm RJ} / T_{\\rm CMB} = \\frac{c^2}{2 \\nu^2 k_B}B^\\prime_\\nu(T_{\\rm CMB, 0}) Here we calculate the RHS of this equation and compare it to the ratio of T_RJ and the result of its transformation to T_CMB. From d424dfdf781df48a52814ad452f9a9bde9c72255 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:29:10 -0700 Subject: [PATCH 13/44] fix: filter out benign numpy warnings https://stackoverflow.com/questions/40845304/runtimewarning-numpy-dtype-size-changed-may-indicate-binary-incompatibility --- pysm3/__init__.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index 0dc674f7..afd79d03 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -1,5 +1,9 @@ # flake8: noqa -from ._astropy_init import * # noqa +from ._astropy_init import * # noqa + +import warnings + +warnings.filterwarnings("ignore", message="numpy.ufunc size changed") try: from importlib.metadata import version, PackageNotFoundError # type: ignore From 2e7e916a560aee597ea6de140022f64af17ecc13 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:34:16 -0700 Subject: [PATCH 14/44] fix codestyle --- pysm3/utils/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pysm3/utils/__init__.py b/pysm3/utils/__init__.py index dee80edb..64ab8dd6 100644 --- a/pysm3/utils/__init__.py +++ b/pysm3/utils/__init__.py @@ -10,8 +10,10 @@ from .data import RemoteData # noqa: F401 import logging + log = logging.getLogger("pysm3") + def has_polarization(m): """Checks if a map or a group of map is polarized @@ -159,7 +161,7 @@ def trapz_step_inplace(freqs, weights, i, m, output): def check_freq_input(freqs): - """ Function to check that the input to `Model.get_emission` is a + """Function to check that the input to `Model.get_emission` is a np.ndarray. This function will convert input scalar frequencies From 3053e5f66e74dc6394e6ad810de2a5d198857cc6 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:42:16 -0700 Subject: [PATCH 15/44] Filter out AstropyDeprecationWarning --- pysm3/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index afd79d03..f6661b41 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -2,7 +2,9 @@ from ._astropy_init import * # noqa import warnings +from astropy.utils.exceptions import AstropyDeprecationWarning +warnings.simplefilter("ignore", category=AstropyDeprecationWarning) warnings.filterwarnings("ignore", message="numpy.ufunc size changed") try: From 79f45554fa5b277c149fa1805133a7cfe8bb2cad Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 21:45:42 -0700 Subject: [PATCH 16/44] filter out all numpy benign warnings --- pysm3/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index f6661b41..bf7937d6 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -5,7 +5,7 @@ from astropy.utils.exceptions import AstropyDeprecationWarning warnings.simplefilter("ignore", category=AstropyDeprecationWarning) -warnings.filterwarnings("ignore", message="numpy.ufunc size changed") +warnings.filterwarnings("ignore", message="may indicate binary incompatibility") try: from importlib.metadata import version, PackageNotFoundError # type: ignore From 917520b86a4cbefab30308e0cdec03cdad46cfd2 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 6 Oct 2021 22:14:46 -0700 Subject: [PATCH 17/44] fix: numpy warning, use float32 instead of float --- pysm3/models/cmb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysm3/models/cmb.py b/pysm3/models/cmb.py index 9aa044bc..ab9d5f4b 100644 --- a/pysm3/models/cmb.py +++ b/pysm3/models/cmb.py @@ -30,7 +30,7 @@ def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ: # do not use normalize weights because that includes a transformation # to spectral radiance and then back to RJ if weights is None: - weights = np.ones(len(freqs), dtype=np.float) + weights = np.ones(len(freqs), dtype=np.float32) scaling_factor = utils.bandpass_unit_conversion( freqs * u.GHz, weights, output_unit=u.uK_RJ, input_unit=u.uK_CMB From 23e875f3de1b327c8f495b3a7f6a96ec8a049fb0 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 07:26:45 -0700 Subject: [PATCH 18/44] fix: use filterwarnings instead of simplefilter --- pysm3/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index bf7937d6..978c8130 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -4,7 +4,7 @@ import warnings from astropy.utils.exceptions import AstropyDeprecationWarning -warnings.simplefilter("ignore", category=AstropyDeprecationWarning) +warnings.filterwarnings("ignore", category=AstropyDeprecationWarning) warnings.filterwarnings("ignore", message="may indicate binary incompatibility") try: From 8e8f2481fe29616d10d354a3d862d92d6fd71851 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 08:05:35 -0700 Subject: [PATCH 19/44] fix: remove unnecessary import --- pysm3/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index 978c8130..777cc07c 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -17,7 +17,6 @@ except PackageNotFoundError: # pragma: no cover __version__ = "unknown" -import sys from .models import * from .sky import Sky from . import units From 137c49d5d3f1117dc3770802884c84b71b2c5cb5 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 08:24:04 -0700 Subject: [PATCH 20/44] fix: remove verbose which is deprecated in healpy 1.15 --- pysm3/models/cmb.py | 2 +- pysm3/models/template.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/pysm3/models/cmb.py b/pysm3/models/cmb.py index ab9d5f4b..7f93d953 100644 --- a/pysm3/models/cmb.py +++ b/pysm3/models/cmb.py @@ -60,7 +60,7 @@ def simulate_tebp_correlated(cl_tebp_arr, nside, lmax, seed): beam_cut = np.ones(3 * nside) for ac in acmb: hp.almxfl(ac, beam_cut, inplace=True) - cmb = np.array(hp.alm2map(acmb, nside, pol=True, verbose=False)) + cmb = np.array(hp.alm2map(acmb, nside, pol=True)) return cmb, aphi diff --git a/pysm3/models/template.py b/pysm3/models/template.py index 006c511c..d85ebad5 100644 --- a/pysm3/models/template.py +++ b/pysm3/models/template.py @@ -132,15 +132,14 @@ def apply_smoothing_and_coord_transform( input_map, lmax=lmax, use_pixel_weights=True if nside > 16 else False, - verbose=False, ) if fwhm is not None: hp.smoothalm( - alm, fwhm=fwhm.to_value(u.rad), verbose=False, inplace=True, pol=True + alm, fwhm=fwhm.to_value(u.rad), inplace=True, pol=True ) if rot is not None: rot.rotate_alm(alm, inplace=True) - smoothed_map = hp.alm2map(alm, nside=nside, verbose=False, pixwin=False) + smoothed_map = hp.alm2map(alm, nside=nside, pixwin=False) else: assert (rot is None) or ( @@ -220,7 +219,7 @@ def read_map(path, nside, unit=None, field=0, map_dist=None): filename = utils.RemoteData().get(path) if (mpi_comm is not None and mpi_comm.rank == 0) or (mpi_comm is None): - output_map = hp.read_map(filename, field=field, verbose=False, dtype=None) + output_map = hp.read_map(filename, field=field, dtype=None) dtype = output_map.dtype # numba only supports little endian if dtype.byteorder == ">": From 586dded8773580825e62d42f7fc02d3ceff0faf8 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 08:26:05 -0700 Subject: [PATCH 21/44] fix: no need of deprecation warning filter anymore --- pysm3/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pysm3/__init__.py b/pysm3/__init__.py index 777cc07c..8f8e2499 100644 --- a/pysm3/__init__.py +++ b/pysm3/__init__.py @@ -4,7 +4,6 @@ import warnings from astropy.utils.exceptions import AstropyDeprecationWarning -warnings.filterwarnings("ignore", category=AstropyDeprecationWarning) warnings.filterwarnings("ignore", message="may indicate binary incompatibility") try: From ef0aef2046844a0705021d1bd97782d1990df758 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 08:28:32 -0700 Subject: [PATCH 22/44] test: try to filter warnings in pytest --- pyproject.toml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 7e7daea4..910f7509 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,3 +5,8 @@ requires = ["setuptools", "wheel"] build-backend = 'setuptools.build_meta' + +[tool.pytest.ini_options] +filterwarnings = [ + 'ignore:may indicate binary incompatibility:RuntimeWarning', +] From e1acb8183d29cd6648daade69746d1619dca1f2b Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 7 Oct 2021 08:41:02 -0700 Subject: [PATCH 23/44] mention warnings to logging and JOSS review in changelog --- CHANGES.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index f07698d4..b8af1da7 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,6 +2,8 @@ ================== - Moved Github repository from `healpy/pysm` to `galsci/pysm`, under the Panexperiment Galactic Science group organization +- Changes in this release are related to the `JOSS Review `_ +- Turned `UserWarning` into `logging`, `Pull Request 88 `_ 3.3.1 (2021-06-30) ================== From c57c3e5b794f24a62a2fdbd8350353ca4db6be33 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 27 Oct 2021 15:05:30 -0700 Subject: [PATCH 24/44] [doc] how to install test requirements From: https://github.com/openjournals/joss-reviews/issues/3783#issuecomment-950407911 --- README.rst | 3 ++- docs/index.rst | 4 ++++ setup.cfg | 3 ++- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/README.rst b/README.rst index 5c6df258..75e9eba5 100644 --- a/README.rst +++ b/README.rst @@ -26,7 +26,8 @@ Install See the `documentation `_ -* Install with ``pip install .`` +* Install with ``pip install .`` or with ``pip install .[test]`` to also install the requirements for running tests +* Optionally, if you have an MPI environment available and you would like to test the MPI capabilities of PySM, install ``mpi4py`` and ``libsharp``, check the documentation link above for more details. * Check code style with ``tox -e codestyle`` * Test with ``pytest`` or ``tox -e test`` * Build docs locally with ``tox -e build_docs`` diff --git a/docs/index.rst b/docs/index.rst index f89766f5..15238669 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -98,6 +98,10 @@ Create a development installation with:: pip install -e . +Install the requirements for testing with:: + + pip install -e .[test] + Execute the unit tests with:: pytest diff --git a/setup.cfg b/setup.cfg index aba0981f..34919a69 100644 --- a/setup.cfg +++ b/setup.cfg @@ -26,7 +26,8 @@ install_requires = [options.extras_require] test = pytest-astropy - mpi4py + pytest + tox docs = sphinx-astropy nbsphinx From b2e3d4b44d3e7bab96f344c92f80a0849436af4c Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 27 Oct 2021 15:08:29 -0700 Subject: [PATCH 25/44] [fix] escape \e to avoid warning --- pysm3/distribution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysm3/distribution.py b/pysm3/distribution.py index 086d91a4..f9a2bb7f 100644 --- a/pysm3/distribution.py +++ b/pysm3/distribution.py @@ -10,7 +10,7 @@ def __init__( In a serial environment, this is only useful if you want to generate a partial sky, pass an array `pixel_indices` with the indices in RING ordering. - in a MPI environment, pass a `mpi4py` communicator and the desired :math:`\ell_{max}` + in a MPI environment, pass a `mpi4py` communicator and the desired :math:`\\ell_{max}` for smoothing and this class will create a ring-based distribution suitable for smoothing with `libsharp`. @@ -23,7 +23,7 @@ def __init__( nside: int Resolution parameter at which this model is to be calculated. smoothing_lmax : int - :math:`\ell_{max}` for the smoothing step, by default :math:`3*N_{side}-1` + :math:`\\ell_{max}` for the smoothing step, by default :math:`3*N_{side}-1` """ self.pixel_indices = pixel_indices self.mpi_comm = mpi_comm From 9bb62280528bad1972d4eea1a9f7dfd44782173c Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 27 Oct 2021 16:18:03 -0700 Subject: [PATCH 26/44] [docs] mention to install pandoc --- README.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/README.rst b/README.rst index 75e9eba5..19cdf411 100644 --- a/README.rst +++ b/README.rst @@ -30,6 +30,7 @@ See the `documentation `_ * Optionally, if you have an MPI environment available and you would like to test the MPI capabilities of PySM, install ``mpi4py`` and ``libsharp``, check the documentation link above for more details. * Check code style with ``tox -e codestyle`` * Test with ``pytest`` or ``tox -e test`` +* Building docs requires ``pandoc``, not the python package, the actual ``pandoc`` command line tool, install it with conda or your package manager * Build docs locally with ``tox -e build_docs`` Support From c2d470feccf406b927d5e121cedbad2eb4e2162a Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 27 Oct 2021 16:25:09 -0700 Subject: [PATCH 27/44] [rtd] install test requirements --- .readthedocs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index 0349ae6f..abcbb6ca 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -11,6 +11,6 @@ python: path: . extra_requirements: - docs - - all + - test formats: [] From 4847b85aec3efff0aa08dd35d9197eb0a104c32b Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Wed, 27 Oct 2021 16:25:33 -0700 Subject: [PATCH 28/44] [rtd] switch to python 3.8 --- .readthedocs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index abcbb6ca..fb063ac0 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -4,7 +4,7 @@ build: image: latest python: - version: 3.7 + version: 3.8 system_packages: true install: - method: pip From 0815cba9740507eb7b9a74987a551e2541eb9145 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 28 Oct 2021 11:36:47 -0700 Subject: [PATCH 29/44] [fix] replace np.int with int --- pysm3/mpi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysm3/mpi.py b/pysm3/mpi.py index ed4adfeb..0a9f9bed 100644 --- a/pysm3/mpi.py +++ b/pysm3/mpi.py @@ -93,10 +93,10 @@ def distribute_rings_libsharp(mpi_comm, nside, lmax): libsharp_grid = libsharp.healpix_grid(nside, rings=local_ring_indices) # returns start index of the ring and number of pixels - startpix, ringpix, _, _, _ = hp.ringinfo(nside, local_ring_indices.astype(np.int64)) + startpix, ringpix, _, _, _ = hp.ringinfo(nside, local_ring_indices.astype(int)) local_npix = libsharp_grid.local_size() - local_pixels = expand_pix(startpix, ringpix, local_npix).astype(np.int) + local_pixels = expand_pix(startpix, ringpix, local_npix).astype(int) local_m_indices = np.arange(mpi_comm.rank, lmax + 1, mpi_comm.size, dtype=np.int32) libsharp_order = libsharp.packed_real_order(lmax, ms=local_m_indices) From c66b16bd8fa64074dedb97b2ce5fc10c49746ff6 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 28 Oct 2021 13:20:52 -0700 Subject: [PATCH 30/44] [docs] better docs for MPI --- docs/mpi.ipynb | 69 +++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 66 insertions(+), 3 deletions(-) diff --git a/docs/mpi.ipynb b/docs/mpi.ipynb index 51ecf8d5..2fcb469c 100644 --- a/docs/mpi.ipynb +++ b/docs/mpi.ipynb @@ -15,14 +15,77 @@ "\n", "The requirements to run with MPI are `mpi4py` and, just for distributed smoothing, `libsharp`.\n", "\n", - "The input maps are read from the first process in order to prevent overloading the filesystem, and then maps are distributed rings by rings across all processes, not overlapping.\n", + "The input maps are read from the first process in order to prevent overloading the filesystem, and then maps are distributed across all processes, not overlapping.\n", "\n", - "The `pysm.MapDistribution` object takes care of handling the metadata about how the data are distributed.\n", + "## Distribution of maps across MPI processes\n", + "\n", + "The `pysm.MapDistribution` object takes care of handling the metadata about how the data are distributed, there are 3 ways to distribute the pixels:\n", + "\n", + "* the most common, which is required to support smoothing over MPI with `libsharp`, is a distribution by ring, where HEALPix rings (stripes of pixels at the same latitude) are distributed symmetrically, i.e. the process that takes first ring close to the North pole also takes the ring close to the South Pole. This is the default distribution if `pixel_indices` is None and a MPI communicator is provided, given that the `libsharp` package is importable:\n", + "\n", + "```python\n", + "map_dist = pysm.MapDistribution(\n", + "pixel_indices=None, nside=nside, mpi_comm=MPI.COMM_WORLD\n", + ")\n", + "```\n", + "\n", + "* in case `libsharp` is not available, it is not possible to smooth the maps, but generating models and integrating bandpass works. The configuration of `MapDistribution` is the same as before, when `libsharp` is not importable. In this case the pixels are distributed uniformly across the processes.\n", + "\n", + "* the final case is a custom distribution, where each MPI process specifies an array of not-overlapping local pixel indices. This is also useful when running serially, if you only want to simulate a small patch of sky:\n", + "\n", + "```python\n", + "map_dist = pysm.MapDistribution(\n", + "pixel_indices=local_pixels_indices, nside=nside, mpi_comm=MPI.COMM_WORLD\n", + ")\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Creation of the `Sky` object\n", + "\n", + "Next you pass the `MapDistribution` object to `Sky` to have all the models being distributed,\n", + "\n", + "```python\n", + "sky = pysm.Sky(nside=nside, preset_strings=[\"d1\", \"s1\", \"f1\", \"a1\"], map_dist=map_dist)\n", + "```\n", "\n", "When the emission is requested, each process independently computes the emission just on its own portion of the sky.\n", "\n", + "```python\n", + "\n", + "m = sky.get_emission(\n", + " freq=np.arange(50, 55) * u.GHz, weights=np.array([0.1, 0.3, 0.5, 0.3, 0.1])\n", + ")[0]\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Distributed smoothing\n", + "\n", "When smoothing is applied, `libsharp` is used to efficiently perform a distributed spherical harmonics transform.\n", "\n", + "```python\n", + "m_smoothed = pysm.apply_smoothing_and_coord_transform(\n", + " m, fwhm=1 * u.deg, map_dist=map_dist\n", + ")\n", + "```\n", + "\n", + "After smoothing, each process holds the smoothed version of their local pixels,\n", + "this can then be used for example to generate timelines, or be collected to a single process for writing (not implemented in PySM currently, please open an issue if you need it)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example\n", + "\n", "Check [`pysm3_mpi.py`](https://github.com/galsci/pysm/blob/master/mpi_examples/pysm3_mpi.py) from the repository as an example.\n", "\n", "Execute with (remove `#` to actually execute it):" @@ -30,7 +93,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ From f55b97e663b4af87858e0f621fa67cc7050aa367 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 28 Oct 2021 13:26:39 -0700 Subject: [PATCH 31/44] [docs] link to MPI tutorial --- docs/index.rst | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/index.rst b/docs/index.rst index 15238669..c5e6c20b 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -50,7 +50,14 @@ Python, which is slower, but easier to debug:: export NUMBA_DISABLE_JIT=1 -In order to run in parallel with MPI, it also needs: +Run PySM with MPI support +------------------------- + +PySM 3 is capable of running across multiple nodes using a MPI communicator. + +See the details in the `MPI section of the tutorial `_. + +In order to run in parallel with MPI, it also needs a functioning MPI environment and: * `mpi4py` From 9ace3e6823948757e227cc721dfef4ad8d1f93c6 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Thu, 28 Oct 2021 13:28:01 -0700 Subject: [PATCH 32/44] [fix] setup mpi requirement --- setup.cfg | 2 ++ tox.ini | 1 + 2 files changed, 3 insertions(+) diff --git a/setup.cfg b/setup.cfg index 34919a69..b7d2f913 100644 --- a/setup.cfg +++ b/setup.cfg @@ -28,6 +28,8 @@ test = pytest-astropy pytest tox +mpi = + mpi4py docs = sphinx-astropy nbsphinx diff --git a/tox.ini b/tox.ini index 03e0ec2e..6a710c20 100644 --- a/tox.ini +++ b/tox.ini @@ -63,6 +63,7 @@ deps = # The following indicates which extras_require from setup.cfg will be installed extras = test + mpi alldeps: all commands = From 6345750e829414d6e2fc989ef816215e2933d459 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Fri, 29 Oct 2021 11:42:30 -0700 Subject: [PATCH 33/44] [fix] fix imports in so pysm model example --- mpi_examples/pysm3_mpi_so_pysm_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mpi_examples/pysm3_mpi_so_pysm_models.py b/mpi_examples/pysm3_mpi_so_pysm_models.py index e633e4ff..0f12c44f 100644 --- a/mpi_examples/pysm3_mpi_so_pysm_models.py +++ b/mpi_examples/pysm3_mpi_so_pysm_models.py @@ -1,7 +1,7 @@ import numpy as np import healpy as hp -import pysm -import pysm.units as u +import pysm3 as pysm +import pysm3.units as u from so_pysm_models import get_so_models import so_pysm_models From c4a16270f94c71bba0c6a91ad957fde1b0846be1 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Fri, 29 Oct 2021 11:48:22 -0700 Subject: [PATCH 34/44] [docs] improve CMBMap --- pysm3/models/cmb.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/pysm3/models/cmb.py b/pysm3/models/cmb.py index 7f93d953..bdbf2db0 100644 --- a/pysm3/models/cmb.py +++ b/pysm3/models/cmb.py @@ -8,9 +8,23 @@ class CMBMap(Model): + """Load one or a set of 3 CMB maps""" + def __init__( self, nside, map_IQU=None, map_I=None, map_Q=None, map_U=None, map_dist=None ): + """ + The input is assumed to be in `uK_CMB` + + Parameters + ---------- + nside: int + HEALPix N_side parameter of the input maps + map_IQU: `pathlib.Path` object + Path to a single IQU map + map_I, map_Q, map_U: `pathlib.Path` object + Paths to the maps to be used as I, Q, U templates. + """ super().__init__(nside=nside, map_dist=map_dist) if map_IQU is not None: self.map = self.read_map(map_IQU, unit=u.uK_CMB, field=(0, 1, 2)) From 5c480eed26946fce164d536b52e1813abeeb7daf Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Fri, 29 Oct 2021 13:29:27 -0700 Subject: [PATCH 35/44] [docs] reference to docs improvement PR --- CHANGES.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGES.rst b/CHANGES.rst index b8af1da7..b823117d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,7 @@ 3.3.2 (unreleased) ================== +- `Improvements to install documentation `_ - Moved Github repository from `healpy/pysm` to `galsci/pysm`, under the Panexperiment Galactic Science group organization - Changes in this release are related to the `JOSS Review `_ - Turned `UserWarning` into `logging`, `Pull Request 88 `_ From 87d8f143a295597c038523cf2c4feebc7a75db3a Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Fri, 29 Oct 2021 13:34:33 -0700 Subject: [PATCH 36/44] [paper] mention tolerance in unit tests --- paper/paper.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/paper/paper.md b/paper/paper.md index db439f3c..e30d1aed 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -77,6 +77,8 @@ The following tables shows the walltime and peak memory usage of this simulation | 1024 | 3m30s 2.3 GB | 7m20s 5.5 GB | | 2048 | 16m10s 8.5 GB | Out of memory | +The models at $N_{side}=512$ have been tested to be equal given a relative tolerance of `1e-5`. + At the moment it is not very useful to run at resolutions higher than $N_{side}=512$ because there is no actual template signal at smaller scales. However, it demonstrates the performance improvements that will make working with higher resolution templates possible. # Future work From 72c26c69649e1877194fe16a67443a50d33f7d20 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Fri, 29 Oct 2021 13:39:27 -0700 Subject: [PATCH 37/44] set release date --- CHANGES.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index b823117d..6218b43d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,4 +1,7 @@ -3.3.2 (unreleased) +3.3.3 (unreleased) +================== + +3.3.2 (2021-10-29) ================== - `Improvements to install documentation `_ From 845fc0d9951db4757ace1c9180b3966a6e2b12be Mon Sep 17 00:00:00 2001 From: "Daniel S. Katz" Date: Mon, 1 Nov 2021 13:30:31 -0500 Subject: [PATCH 38/44] small changes in paper --- paper/paper.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/paper/paper.md b/paper/paper.md index e30d1aed..f40d80cc 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -37,25 +37,25 @@ bibliography: paper.bib # Statement of Need The Cosmic Microwave Background (CMB) radiation, emitted just 370 thousand years after the Big Bang, is a pristine probe of the Early Universe. After being emitted at high temperatures, the CMB was redshifted by the subsequent 13.8 billion years of cosmic expansion, such that it is brightest at microwave frequencies today. -However, our own Milky Way galaxy also emits in the microwave portion of the spectrum, obscuring our view of the CMB. Examples of this emission are thermal radiation by interstellar dust grains, and synchrotron emission by relativistic electrons spiraling in magnetic fields. +However, our own Milky Way galaxy also emits in the microwave portion of the spectrum, obscuring our view of the CMB. Examples of this emission are thermal radiation by interstellar dust grains and synchrotron emission by relativistic electrons spiraling in magnetic fields. Cosmologists need to create synthetic maps of the CMB and of the galactic emission based on available data and on physical models that extrapolate observations to different frequencies. The resulting maps are useful to test data reduction algorithms, to understand residual systematics, to forecast maps produced by future instruments, to run Monte Carlo analysis for noise estimation, and more. # Summary -The Python Sky Model (PySM) is a Python package used by Cosmic Microwave Background (CMB) experiments to simulate maps, in HEALPix [@gorski05; @healpy09] pixelization, of the various diffuse astrophysical components of Galactic emission relevant at CMB frequencies (i.e. dust, synchrotron, free-free and Anomalous Microwave Emission), as well as the CMB itself. These maps may be integrated over a given instrument bandpass and smoothed with a given instrument beam. +The Python Sky Model (PySM) is a Python package used by Cosmic Microwave Background (CMB) experiments to simulate maps, in HEALPix [@gorski05; @healpy09] pixelization, of the various diffuse astrophysical components of Galactic emission relevant at CMB frequencies (i.e., dust, synchrotron, free-free and Anomalous Microwave Emission), as well as the CMB itself. These maps may be integrated over a given instrument bandpass and smoothed with a given instrument beam. The template emission maps used by PySM are based on Planck [@planck18] and WMAP [@wmap13] data and are noise-dominated at small scales. Therefore, PySM simulation templates are smoothed to retain the large-scale information, and then supplemented with modulated Gaussian realizations at smaller scales. This strategy allows one to simulate data at higher resolution than the input maps. -PySM 2 [@pysm17], released in 2016, has become the de-facto standard for simulating Galactic emission, for example it is used by CMB-S4, Simons Observatory, LiteBird, PICO, CLASS, POLARBEAR and other CMB experiments, as shown by the [80+ citations of the PySM 2 publication](https://scholar.google.com/scholar?start=0&hl=en&as_sdt=2005&sciodt=0,5&cites=16628417670342266167&scipsc=). +PySM 2 [@pysm17], released in 2016, has become the de-facto standard for simulating Galactic emission; it is used, for example, by CMB-S4, Simons Observatory, LiteBird, PICO, CLASS, POLARBEAR, and other CMB experiments, as shown by the [80+ citations of the PySM 2 publication](https://scholar.google.com/scholar?start=0&hl=en&as_sdt=2005&sciodt=0,5&cites=16628417670342266167&scipsc=). As the resolution of upcoming experiments increases, the PySM 2 software has started to show some limitations: * Emission templates are provided at 7.9 arcminutes resolution (HEALPix $N_{side}=512$), while the next generation of CMB experiments will require sub-arcminute resolution. -* The software is implemented in pure `numpy`, meaning that it has significant memory overhead and is not multi-threaded, precluding simply replacing the current templates with higher-resolution versions -* Emission templates are included in the PySM 2 Python package, this is still practical when each of the roughly 40 input maps is ~10 Megabytes, but will not be if they are over 1 Gigabyte. +* The software is implemented in pure `numpy`, meaning that it has significant memory overhead and is not multi-threaded, precluding simply replacing the current templates with higher-resolution versions. +* Emission templates are included in the PySM 2 Python package, which is still practical when each of the roughly 40 input maps is ~10 Megabytes, but will not be if they are over 1 Gigabyte. The solution to these issues was to reimplement PySM from scratch focusing of these features: * Reimplement all the models with the `numba` [@numba] Just-In-Time compiler for Python to reduce memory overhead and optimize performance: the whole integration loop of a template map over the frequency response of an instrument is performed in a single pass in automatically compiled and multi-threaded Python code. -* Use MPI through `mpi4py` to coordinate execution of PySM 3 across multiple nodes, this allows to support template maps at a resolution up to 0.4 arcminutes (HEALPix $N_{side}=8192$). +* Use MPI through `mpi4py` to coordinate execution of PySM 3 across multiple nodes, this allows supporting template maps at a resolution up to 0.4 arcminutes (HEALPix $N_{side}=8192$). * Rely on `libsharp` [@libsharp], a distributed implementation of spherical harmonic transforms, to smooth the maps with the instrument beam when maps are distributed over multiple nodes with MPI. * Employ the data utilities infrastructure provided by `astropy` [@astropy2013; @astropy2018] to download the input templates and cache them when requested. @@ -79,11 +79,11 @@ The following tables shows the walltime and peak memory usage of this simulation The models at $N_{side}=512$ have been tested to be equal given a relative tolerance of `1e-5`. -At the moment it is not very useful to run at resolutions higher than $N_{side}=512$ because there is no actual template signal at smaller scales. However, it demonstrates the performance improvements that will make working with higher resolution templates possible. +At the moment it is not very useful to run at resolutions higher than $N_{side}=512$ because there is no actual template signal at smaller scales. However, this demonstrates the performance improvements that will make working with higher resolution templates possible. # Future work -PySM 3 opens the way to implement a new category of models at much higher resolution. However, instead of just upgrading the current models to smaller scales we want to also update them with the latest knowledge of Galactic emission and gather feedback from each of the numerous CMB experiments. For this reason we are collaborating with the Panexperiment Galactic Science group to lead the development of the new class of models to be included in PySM 3. +PySM 3 opens the way to implement a new category of models at much higher resolution. However, instead of just upgrading the current models to smaller scales, we want to also update them with the latest knowledge of Galactic emission and gather feedback from each of the numerous CMB experiments. For this reason we are collaborating with the Panexperiment Galactic Science group to lead the development of the new class of models to be included in PySM 3. # How to cite From b07a27b5c37f87e72aa001f32bc47b0cdd03aeed Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 1 Nov 2021 11:47:37 -0700 Subject: [PATCH 39/44] [paper] replace journal shortcuts --- paper/paper.bib | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/paper/paper.bib b/paper/paper.bib index f1c845ff..0fa6a2db 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -3,7 +3,7 @@ @ARTICLE{gorski05 {Wandelt}, B.~D. and {Hansen}, F.~K. and {Reinecke}, M. and {Bartelmann}, M.}, title = "{HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere}", - journal = {\apj}, + journal = {The Astrophysical Journal}, eprint = {astro-ph/0409513}, keywords = {Cosmology: Cosmic Microwave Background, Cosmology: Observations, Methods: Statistical}, year = 2005, @@ -51,7 +51,7 @@ @ARTICLE{wmap13 {Kogut}, A. and {Limon}, M. and {Meyer}, S.~S. and {Tucker}, G.~S. and {Wright}, E.~L.}, title = "{Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Final Maps and Results}", - journal = {\apjs}, + journal = {The Astrophysical Journal Supplement Series}, archivePrefix = "arXiv", eprint = {1212.5225}, keywords = {cosmic background radiation, cosmology: observations, dark matter, early universe, instrumentation: detectors, space vehicles, space vehicles: instruments, telescopes}, @@ -107,7 +107,7 @@ @article{libsharp ISSN={1432-0746}, url={http://dx.doi.org/10.1051/0004-6361/201321494}, DOI={10.1051/0004-6361/201321494}, - journal={Astronomy & Astrophysics}, + journal={Astronomy and Astrophysics}, publisher={EDP Sciences}, author={Reinecke, M. and Seljebotn, D. S.}, year={2013}, @@ -122,7 +122,7 @@ @article{astropy2013 Doi = {10.1051/0004-6361/201322068}, Eid = {A33}, Eprint = {1307.6212}, -Journal = {\aap}, +Journal = {Astronomy and Astrophysics}, Keywords = {methods: data analysis, methods: miscellaneous, virtual observatory tools}, Month = oct, Pages = {A33}, @@ -177,7 +177,7 @@ @ARTICLE{astropy2018 {Weaver}, B.~A. and {Whitmore}, J.~B. and {Woillez}, J. and {Zabalza}, V. and {Astropy Contributors}}, title = "{The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package}", - journal = {\aj}, + journal = {The Astronomical Journal}, keywords = {methods: data analysis, methods: miscellaneous, methods: statistical, reference systems, Astrophysics - Instrumentation and Methods for Astrophysics}, year = 2018, month = sep, From 61363ab7b9c71fe5160994e68abea2acbe65dbaf Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 1 Nov 2021 12:05:21 -0700 Subject: [PATCH 40/44] [paper] fix bibliography capitalization --- paper/paper.bib | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/paper/paper.bib b/paper/paper.bib index 0fa6a2db..5756590d 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -24,12 +24,12 @@ @article{healpy09 number = {35}, pages = {1298}, author = {Andrea Zonca and Leo P. Singer and Daniel Lenz and Martin Reinecke and Cyrille Rosset and Eric Hivon and Krzysztof M. Gorski}, - title = {healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python}, + title = {{healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python}}, journal = {Journal of Open Source Software} } @article{pysm17, - title={The Python Sky Model: software for simulating the Galactic microwave sky}, + title={{The Python Sky Model: software for simulating the Galactic microwave sky}}, volume={469}, ISSN={1365-2966}, url={http://dx.doi.org/10.1093/mnras/stx949}, From cec3de6e9bee8d24d9369f6f3f2bae36610e1a89 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 1 Nov 2021 14:29:51 -0700 Subject: [PATCH 41/44] [paper] update CITATION to cite JOSS instead of Arxiv --- CITATION | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/CITATION b/CITATION index 0fc4176f..74f46404 100644 --- a/CITATION +++ b/CITATION @@ -1,13 +1,18 @@ -@misc{pysm3, - title={The Python Sky Model 3 software}, - author={Andrea Zonca and Ben Thorne and Nicoletta Krachmalnicoff and Julian Borrill}, - year={2021}, - eprint={2108.01444}, - archivePrefix={arXiv}, - primaryClass={astro-ph.IM} +@article{Zonca_2021, + doi = {10.21105/joss.03783}, + url = {https://doi.org/10.21105/joss.03783}, + year = {2021}, + publisher = {The Open Journal}, + volume = {6}, + number = {67}, + pages = {3783}, + author = {Andrea Zonca and Ben Thorne and Nicoletta Krachmalnicoff and Julian Borrill}, + title = {{The Python Sky Model 3 software}}, + journal = {Journal of Open Source Software} } + @article{Thorne_2017, - title={The Python Sky Model: software for simulating the Galactic microwave sky}, + title={{The Python Sky Model: software for simulating the Galactic microwave sky}}, volume={469}, ISSN={1365-2966}, url={http://dx.doi.org/10.1093/mnras/stx949}, From d300ba9ed6e8225ae8612c09bc3a9884948467ea Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 31 Jan 2022 22:33:04 -0800 Subject: [PATCH 42/44] test: implement unit test for remotedata --- pysm3/tests/test_utils.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/pysm3/tests/test_utils.py b/pysm3/tests/test_utils.py index 0cd920bc..4a7d11c7 100644 --- a/pysm3/tests/test_utils.py +++ b/pysm3/tests/test_utils.py @@ -31,3 +31,13 @@ def test_bandpass_unit_conversion(): ) for pol in [0, 1]: assert_quantity_allclose(expected_map[pol], CMB_thermo_int[pol], rtol=1e-4) + +def test_remotedata(tmp_path): + import os + data_folder = tmp_path / "data" + data_folder.mkdir() + test_file = data_folder / "testfile.txt" + test_file.touch() + os.environ["PYSM_LOCAL_DATA"] = str(data_folder) + filename = pysm3.utils.RemoteData().get("testfile.txt") + assert filename == str(test_file) From a86cc73b4fa23e07937bfc2ba9f7c20c6a64a47c Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 31 Jan 2022 22:43:08 -0800 Subject: [PATCH 43/44] test: enable logging in pytest --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 910f7509..5fb4b5c7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,3 +10,5 @@ build-backend = 'setuptools.build_meta' filterwarnings = [ 'ignore:may indicate binary incompatibility:RuntimeWarning', ] +log_cli = true +log_cli_level = "DEBUG" From b30f28fd091b5dfb3c603ace18a00ab1d2618e35 Mon Sep 17 00:00:00 2001 From: Andrea Zonca Date: Mon, 31 Jan 2022 22:45:58 -0800 Subject: [PATCH 44/44] run test with global path --- pysm3/tests/test_utils.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pysm3/tests/test_utils.py b/pysm3/tests/test_utils.py index 4a7d11c7..21fe20d1 100644 --- a/pysm3/tests/test_utils.py +++ b/pysm3/tests/test_utils.py @@ -41,3 +41,12 @@ def test_remotedata(tmp_path): os.environ["PYSM_LOCAL_DATA"] = str(data_folder) filename = pysm3.utils.RemoteData().get("testfile.txt") assert filename == str(test_file) + + +def test_remotedata_globalpath(tmp_path): + data_folder = tmp_path / "data" + data_folder.mkdir() + test_file = data_folder / "testfile.txt" + test_file.touch() + filename = pysm3.utils.RemoteData().get(str(test_file)) + assert filename == str(test_file)