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?
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 0x7fbf5cb43150>

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 0x7fbf5d0ed910>
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")

Exercise

Create a PySMSky object with pure CMB (model c1), in uK_CMB

Get the CMB signal for 2 channels at 2 different frequencies and verify that the maps are the same to double precision error (~1e-14). To specify a channel at a single frequency instead of a bandpass, you can specify:

bandpasses={"ch0":(100,1)}

Where the first element of the tuple is frequency in GHz, the second its weight.

In [ ]: