Simulated Sky Signal in map domain

This lesson is about simulating the input sky signal using PySM 3.

PySM 3

If you used PySM in the past, you most probably used PySM 2 from https://github.com/bthorne93/PySM_public.

PySM 3 is a rewrite of PySM which offers all the same functionality and the same models of PySM 2 but is focused on:

  • improving performance using just-in-time compilation and multi-threading with numba
  • lowering memory requirements by reworking the underlying algorithms
  • improved capability of running in parallel with MPI

It is available from https://github.com/healpy/pysm, it is still missing a few models and the documentation but is already integrated into TOAST to overcame the strong performance limits of PySM 2.

If anyone is interested in learning more about PySM 3, check the PySM 3 tutorial, we can work through this during the hack day.

PySMSky

The lower level TOAST class is PySMSky, it performs the following operations:

  • initialize PySM with the input sky configuration
  • loop through all channels and for each calls PySM to generate the sky emission at all frequencies in the bandpass and integrate
In [1]:
# Load common tools for all lessons
import sys
sys.path.insert(0, "..")
from lesson_tools import (
    fake_focalplane
)

# Capture C++ output in the jupyter cells
%reload_ext wurlitzer
In [2]:
import toast
import healpy as hp
import numpy as np
In [3]:
env = toast.Environment.get()                         
env.set_log_level("DEBUG")
In [4]:
focal_plane = fake_focalplane()
In [5]:
from toast.todmap import PySMSky
PySMSky?

Init signature:
PySMSky(
    comm=None,
    pixels='pixels',
    out='pysm',
    nside=None,
    pysm_sky_config=None,
    pysm_precomputed_cmb_K_CMB=None,
    pysm_component_objects=None,
    init_sky=True,
    pixel_indices=None,
    units='K_CMB',
    map_dist=None,
)
Docstring:     
Create a bandpass-integrated sky map with PySM

Requires PySM 3. It initializes the `pysm.Sky` object either
in the constructor (`init_sky=True`) or when the `exec` method
is executed (`init_sky=False`).
Inizialization of the sky will load templates from the first process
of the MPI communicator, copy to all processes and then select the local
rings, distributed as required by `libsharp` for eventual smoothing later.
If another pixel distribution is required, it can be specified with `pixel_indices`,
however a different pixel distribution can only be used if performing no smoothing.

Args:
    comm (mpi4py.MPI.Comm): MPI communicator
    pixels (str): the name of the cache object (<pixels>_<detector>)
        containing the pixel indices to use.
    out (str): accumulate data to the cache with name <out>_<detector>.
        If the named cache objects do not exist, then they are created.
    nside (int): :math:`N_{side}` for PySM
    pysm_sky_config (list(str)): list of PySM components, e.g. ["s1", "d1"],
        this will be passed as `preset_strings` to `pysm.Sky`
    pysm_precomputed_cmb_K_CMB (str): obsolete, will be removed shortly
    pysm_component_objects (list(pysm.Model)): extra sky components that
        inherits from `pysm.Model` to be passed to `pysm.Sky` as `components_objects`
    init_sky (bool): Initializes the sky in the constructor if True, in `exec` if False
    pixel_indices (np.array(int)): List of pixel indices, use None to use the standard
        ring-based libsharp distribution
    units (str): Output units.
File:           /global/common/software/cmb/cori/cmbenv-intel_20191012/cmbenv_aux/lib/python3.6/site-packages/toast/todmap/pysm.py
Type:           type
Subclasses:
In [6]:
NSIDE = 64
npix = hp.nside2npix(NSIDE)

PySM models

You can find out details about all the models available in PySM 2 and 3 at:

https://pysm-public.readthedocs.io/en/latest/models.html

In [7]:
pysm_sky_config = ["s1", "f1", "a1", "d1"]
In [8]:
pysm_sky = PySMSky(comm=None,
       pixel_indices=None,
       nside=NSIDE,
       units="uK_RJ",
       pysm_sky_config=pysm_sky_config
)
In [9]:
pysm_sky
Out[9]:
<toast.todmap.pysm.PySMSky at 0x2aaaeae8c5f8>

PySM Sky object

We can directly access the underlying PySM.Sky object as the sky attribute of the PySMSky object.

In [10]:
pysm_sky.sky
Out[10]:
<pysm.sky.Sky at 0x2aaaeae8c7f0>
In [11]:
import pysm.units as u
In [12]:
pysm_sky.sky.get_emission(12 * u.GHz)
Out[12]:
$[[594.75165,~477.8017,~464.08673,~\dots,~423.69739,~490.87775,~478.81104],~ [-10.309416,~-20.71452,~-28.465218,~\dots,~-1.3051515,~-13.430266,~0.012164834],~ [-14.353321,~15.995528,~4.725172,~\dots,~-20.496,~11.566074,~-23.161991]] \; \mathrm{\mu K_{{RJ}}}$
In [13]:
%matplotlib inline
## _ refers to the output of the previous cell, so this works only if you run cells in sequence
hp.mollview(_[0], cmap="coolwarm", min=-100, max=1e4)

Execute the PySMSky object

First we need bandpasses for the channels, first element of the tuple is frequency, second element are weights, we define a top-hat of 10 points with a bandwidth of 10 GHz:

In [14]:
bandpasses = {}
for ch in focal_plane: # loops through dict keys
    bandpasses[ch] = (np.linspace(65, 75, 10), np.ones(10))
In [15]:
bandpasses
Out[15]:
{'0A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '0B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '1A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '1B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '2A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '2B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '3A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '3B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '4A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '4B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '5A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '5B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '6A': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])),
 '6B': (array([65.        , 66.11111111, 67.22222222, 68.33333333, 69.44444444,
         70.55555556, 71.66666667, 72.77777778, 73.88888889, 75.        ]),
  array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]))}
In [16]:
local_maps = {}
pysm_sky.exec(local_maps, out="sky", bandpasses=bandpasses)
In [17]:
hp.mollview(local_maps["sky_0A"][0], cmap="coolwarm", min=0, max=1e3, unit="uK_RJ")
In [18]:
local_maps["sky_0B"][0]-local_maps["sky_0A"][0]
Out[18]:
array([0., 0., 0., ..., 0., 0., 0.], dtype=float32)
In [19]:
bandpasses["0B"] = (np.linspace(63, 73, 10), np.ones(10))
In [20]:
local_maps = {}
pysm_sky.exec(local_maps, out="sky", bandpasses=bandpasses)
In [21]:
hp.mollview(local_maps["sky_0A"][0]-local_maps["sky_0B"][0], cmap="coolwarm", unit="uK_RJ")
In [22]:
hp.gnomview(local_maps["sky_0A"][0]-local_maps["sky_0B"][0], rot=(0,0), xsize=5000, ysize=2000, cmap="coolwarm")
In [ ]: