Output car and modify Nside after smoothing

.github/workflows/ci_tests.yml

@@ -22,25 +22,25 @@ jobs:
     strategy:
       matrix:
         include:
-          - name: Code style checks
-            os: ubuntu-latest
-            python: 3.x
-            toxenv: codestyle
+           - name: Code style checks
+             os: ubuntu-latest
+             python: 3.x
+             toxenv: codestyle
 
-          - name: Python 3.7 with minimal dependencies
-            os: ubuntu-latest
-            python: 3.7
-            toxenv: py37-test
+           - name: Python 3.7 with minimal dependencies
+             os: ubuntu-latest
+             python: 3.7
+             toxenv: py37-test
 
-          - name: Python 3.8 with all optional dependencies and coverage checking
-            os: ubuntu-latest
-            python: 3.8
-            toxenv: py38-test-alldeps-cov
+           - name: Python 3.8 with all optional dependencies and coverage checking
+             os: ubuntu-latest
+             python: 3.8
+             toxenv: py38-test-alldeps-cov
 
-          - name: OS X - Python 3.8
-            os: macos-latest
-            python: 3.8
-            toxenv: py38-test
+             # - name: OS X - Python 3.8
+             #   os: macos-10.15
+             #   python: 3.8
+             #   toxenv: py38-test
 
           # - name: Python 3.7 with oldest supported version of all dependencies
           #   os: ubuntu-16.04

CHANGES.rst

@@ -1,8 +1,8 @@
 3.4.0 (unreleased)
 ==================
 
+- `apply_smoothing_and_coord_transform` now supports a different output resolution and supports doing both HEALPix and CAR in the same execution, also added best practices for dealing with resolution in the docs `PR 125 <https://github.com/galsci/pysm/pull/125>`_
 - Model has `max_nside` attribute which specifies its max resolution `PR 124 <https://github.com/galsci/pysm/pull/124>`_
-- Running unit tests (except MPI) on MAC OS X in Github actions `PR 122 <https://github.com/galsci/pysm/pull/122>`_
 
 3.4.0b3 (2022-03-28)
 ====================

docs/index.rst

@@ -28,6 +28,19 @@ For example free-free:
 
   models
 
+Best practices for model execution
+----------------------------------
+
+PySM 3, in order to have the same behaviour of PySM 2, uses `hp.ud_grade` to change the resolution of the map if the HEALPix $N_{side}$ of the input templates is different from the requested output resolution, specified in the emission class (subclass of :py:class:`~pysm3.Model`) or in the :py:class:`~pysm3.Sky` class.
+
+`hp.ud_grade` generally creates artifacts in the spectra and should be avoided unless the specific application can tolerate that.
+
+Therefore we recommend to execute all PySM 2 derived models (e.g. `d0` to `d8`) at their native $N_{side}$ of 512, and then use the `output_nside` parameter of :py:func:`apply_smoothing_and_coord_transform` to transform to the target resolution, whether higher or lower, in Spherical Harmonics domain.
+
+PySM 3 native models (e.g. `d9` to `d11` and `s4` to `s6`) instead have precomputed templates (generated in Spherical Harmonics domain) from $N_{side}$ 2048 to 8192. Therefore we recommend to execute them at 2048 if the target $N_{side}$ is 1024 or lower and at 2*$N_{side}$ if the target is 2048 or higher, except 8192 which is the highest possible resolution.
+
+There are some exceptions, so always check in the documentation or the `max_nside` parameter of the models, for example the `d12` model, even if native to PySM 3, uses externally provided maps which are only available at $N_{side}$=2048.
+
 Dependencies
 ============
 

pysm3/models/spdust.py

@@ -7,7 +7,7 @@
 
 
 class SpDust(Model):
-    """ Implementation of the SpDust2 code of (Ali-Haimoud et al 2012)
+    """Implementation of the SpDust2 code of (Ali-Haimoud et al 2012)
     evaluated for a Cold Neutral Medium.
     See:
     * https://arxiv.org/abs/1003.4732
@@ -26,7 +26,7 @@ def __init__(
         unit_I=None,
         map_dist=None,
     ):
-        """ This function initializes the spinning dust model
+        """This function initializes the spinning dust model
 
         Parameters
         ----------
@@ -149,7 +149,7 @@ def get_emission(self, freqs: u.GHz, weights=None):
                 self.freq_ref_I.value,
                 self.freq_peak.value,
                 self.emissivity,
-                self.pol_angle,
+                self.pol_angle.value,
                 self.pol_frac,
             )
             << u.uK_RJ

pysm3/models/template.py

@@ -16,6 +16,13 @@
 from .. import mpi
 import gc
 
+try:
+    import pixell.enmap
+    import pixell.curvedsky
+    import pixell.sharp
+except ImportError:
+    pixell = None
+
 log = logging.getLogger("pysm3")
 
 
@@ -42,6 +49,9 @@ def __init__(self, nside, max_nside=None, map_dist=None):
             MPI communicator object (optional, default=None).
         nside: int
             Resolution parameter at which this model is to be calculated.
+        max_nside: int
+            Keeps track of the the maximum Nside this model is available at
+            by default 512 like PySM 2 models
         smoothing_lmax : int
             :math:`\\ell_{max}` for the smoothing step, by default :math:`2*N_{side}`
         """
@@ -112,7 +122,16 @@ def get_emission(self, freqs: u.GHz, weights=None) -> u.uK_RJ:
 
 
 def apply_smoothing_and_coord_transform(
-    input_map, fwhm=None, rot=None, lmax=None, map_dist=None
+    input_map,
+    fwhm=None,
+    rot=None,
+    lmax=None,
+    output_nside=None,
+    output_car_resol=None,
+    return_healpix=True,
+    return_car=False,
+    input_alm=False,
+    map_dist=None,
 ):
     """Apply smoothing and coordinate rotation to an input map
 
@@ -125,7 +144,7 @@ def apply_smoothing_and_coord_transform(
     input_map : ndarray
         Input map, of shape `(3, npix)`
         This is assumed to have no beam at this point, as the
-        simulated small scale tempatle on which the simulations are based
+        simulated small scale template on which the simulations are based
         have no beam.
     fwhm : astropy.units.Quantity
         Full width at half-maximum, defining the
@@ -134,35 +153,103 @@ def apply_smoothing_and_coord_transform(
         Apply a coordinate rotation give a healpy `Rotator`, e.g. if the
         inputs are in Galactic, `hp.Rotator(coord=("G", "C"))` rotates
         to Equatorial
+    output_nside : int
+        HEALPix output map Nside, if None, use the same as the input
+    lmax : int
+        lmax for the map2alm step, if None, it is set to 2.5 * nside
+        if output_nside is equal or higher than nside.
+        It is set to 1.5 * nside if output_nside is lower than nside
+    output_car_resol : astropy.Quantity
+        CAR output map resolution, generally in arcmin
+    return_healpix : bool
+        Whether to return the HEALPix map
+    return_car : bool
+        Whether to return the CAR map
+    input_alm : np.array
+        Instead of starting from a map, `input_map` is a set of Alm
 
     Returns
     -------
-    smoothed_map : np.ndarray
-        Array containing the smoothed sky
+    smoothed_map : np.ndarray or tuple of np.ndarray
+        Array containing the smoothed sky or tuple of HEALPix and CAR maps
     """
 
-    if map_dist is None:
+    if not input_alm:
         nside = hp.get_nside(input_map)
-        alm = hp.map2alm(
-            input_map,
-            lmax=lmax,
-            use_pixel_weights=True if nside > 16 else False,
-        )
+        if output_nside is None:
+            output_nside = nside
+
+    if hasattr(input_map, "unit"):
+        unit = input_map.unit
+    else:
+        unit = 1
+
+    if lmax is None:
+        if nside == output_nside:
+            lmax = int(2.5 * output_nside)
+        elif output_nside > nside:
+            lmax = int(2.5 * nside)
+        elif output_nside < nside:
+            lmax = int(1.5 * nside)
+
+    output_maps = []
+
+    if map_dist is None:
+        if input_alm:
+            alm = input_map.copy()
+        else:
+            if lmax <= 1.5 * nside:
+                log.info("Using map2alm with pixel weights")
+                alm = hp.map2alm(
+                    input_map,
+                    lmax=lmax,
+                    use_pixel_weights=True if nside > 16 else False,
+                )
+            else:
+                alm, error, n_iter = hp.map2alm_lsq(
+                    input_map, lmax=lmax, mmax=lmax, tol=1e-6, maxiter=100
+                )
+                log.info(
+                    "Used map2alm_lsq, converged in %d iterations, residual relative error %.2g",
+                    n_iter,
+                    error,
+                )
+                if n_iter == 100:
+                    log.warning(
+                        "hp.map2alm_lsq did not converge in %d iterations,"
+                        + " residual relative error is %.2g",
+                        n_iter,
+                        error,
+                    )
         if fwhm is not None:
             hp.smoothalm(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, pixwin=False)
-
+        if return_healpix:
+            if input_alm:
+                assert (
+                    output_nside is not None
+                ), "If inputting Alms, specify output_nside"
+            output_maps.append(hp.alm2map(alm, nside=output_nside, pixwin=False) * unit)
+        if return_car:
+            shape, wcs = pixell.enmap.fullsky_geometry(
+                output_car_resol.to_value(u.radian), dims=(3,)
+            )
+            ainfo = pixell.sharp.alm_info(lmax=lmax)
+            output_maps.append(
+                pixell.curvedsky.alm2map(
+                    alm, pixell.enmap.empty(shape, wcs), ainfo=ainfo
+                )
+                * unit
+            )
     else:
         assert (rot is None) or (
             rot.coordin == rot.coordout
         ), "No rotation supported in distributed smoothing"
-        smoothed_map = mpi.mpi_smoothing(input_map, fwhm, map_dist)
+        output_maps.append(mpi.mpi_smoothing(input_map, fwhm, map_dist))
+        assert not return_car, "No CAR output supported in Libsharp smoothing"
 
-    if hasattr(input_map, "unit"):
-        smoothed_map <<= input_map.unit
-    return smoothed_map
+    return output_maps[0] if len(output_maps) == 1 else tuple(output_maps)
 
 
 def apply_normalization(freqs, weights):
@@ -343,7 +430,7 @@ def read_txt(path, mpi_comm=None, **kwargs):
 
 
 def read_alm(path, has_polarization=True, unit=None, map_dist=None):
-    """Read :math:`a_{\ell m}` from a FITS file
+    """Read :math:`a_{\\ell m}` from a FITS file
 
     Parameters
     ----------
@@ -352,8 +439,8 @@ def read_alm(path, has_polarization=True, unit=None, map_dist=None):
     has_polarization : bool
         read only temperature alm from file or also polarization
     map_dist : pysm.MapDistribution
-        :math:`\ell_{max}` should be the same of the :math:`\ell_{max}` in the file
-        and :math:`m_{max}=\ell_{max}`.
+        :math:`\\ell_{max}` should be the same of the :math:`\\ell_{max}` in the file
+        and :math:`m_{max}=\\ell_{max}`.
     """
 
     filename = utils.RemoteData().get(path)
@@ -377,7 +464,7 @@ def read_alm(path, has_polarization=True, unit=None, map_dist=None):
 
 
 def read_cl(path, has_polarization=True, unit=None, map_dist=None):
-    """Read :math:`a_{\ell m}` from a FITS file
+    """Read :math:`a_{\\ell m}` from a FITS file
 
     Parameters
     ----------
@@ -386,8 +473,8 @@ def read_cl(path, has_polarization=True, unit=None, map_dist=None):
     has_polarization : bool
         read only temperature alm from file or also polarization
     map_dist : pysm.MapDistribution
-        :math:`\ell_{max}` should be the same of the :math:`\ell_{max}` in the file
-        and :math:`m_{max}=\ell_{max}`.
+        :math:`\\ell_{max}` should be the same of the :math:`\\ell_{max}` in the file
+        and :math:`m_{max}=\\ell_{max}`.
     """
 
     filename = utils.RemoteData().get(path)

pysm3/sky.py

@@ -60,6 +60,9 @@ def create_components_from_config(config, nside, map_dist=None):
 class Sky(Model):
     """Sky is the main interface to PySM
 
+    Please read the 'Best practices for model execution' section in the
+    documentation homepage before running PySM 3 models.
+
     It accepts the configuration of the desired components in 3 different
     ways: `preset_strings`, `component_config` or `component_objects`,
     see details below.
@@ -95,6 +98,9 @@ def __init__(
         nside : int
             Requested output NSIDE, inputs will be degraded
             using :func:`healpy.ud_grade`
+        max_nside: int
+            Keeps track of the the maximum Nside this model is available at
+            by default 512 like PySM 2 models
         preset_strings : list of str
             List of strings identifiers for the models included in PySM 3,
             these are exactly the same models included in PySM 2, e.g.
@@ -146,7 +152,7 @@ def add_component(self, component):
         self.components.append(component)
 
     def get_emission(self, freq, weights=None):
-        """ This function returns the emission at a frequency, set of
+        """This function returns the emission at a frequency, set of
         frequencies, or over a bandpass.
         """
         output = self.components[0].get_emission(freq, weights=weights)

pysm3/tests/test_dust_pysm3.py

@@ -73,7 +73,7 @@ def test_gnilc_857(model_tag):
         )
     )
     scaling = (freq / freq_ref) ** (beta - 2)
-    scaling *= blackbody_ratio(freq, freq_ref, Td.to_value(u.K))
+    scaling *= blackbody_ratio(freq.value, freq_ref.value, Td.to_value(u.K))
 
     assert_quantity_allclose(input_template * scaling, output, rtol=1e-6)