Store templates in Alm instead of maps?

[zonca]

Problem

PySM needs to store the templates at very different resolutions, we want people to be able to do maps at low nside without the need of downloading and processing very high resolution maps.

Current implementation

Currently PySM 3 implements the same method of PySM 2. Input templates are stored at a fixed N_side (512 for PySM 2 models), higher resolution models implemented in so_pysm_models have both 512 and 4096. Any other N_side is computed via ud_grade from the higher resolution map. ud_grade is quick but heavily smoothes the spectrum, in fact both SO and S4 run simulations at 512 and 4096, so effectively never suffering this.

See for example impact of ud_grade on a power law signal with a slope of -2:

"Provide all the maps"

The simplest option is to pre-compute maps using alm2map after clipping the Alms to 3 * N_side - 1:

It is a lot of maps but the N_side 8192 is dominating the storage space anyway, the sum of all other maps is less than half of it. We should provide maps at all N_sides from 8192 down to 128 (Lower resolution for SO)

In this case, we download, cache and read in memory the maps at the requested N_side, then do bandpass integration, finally do beam smoothing. So if we are simulating 1 channel at a time we do 1 map2alm and 1 alm2map. Just as a reference, at 2048 on Jupyter@NERSC map2alm on a IQU map takes 80s, alm2map takes 40s, so about 120s.

"Provide the Alms"

Alternatively we could provide the Alms instead of maps, that would be at variable resolution, for example from 8192 (quoting N_side, I assume ell_max is 3*N_side-1) down to 512. Then we first apply the beam smoothing to all templates with almxfl, then transform to maps with alm2map at the target N_side. So the Sky object now would also have a specific beam instead of having only a specific N_side, but I guess most people run PySM for 1 component and 1 channel at a time. Most components need to multiply maps in pixel space, so we need to transform the templates to maps before running the model.

Then we do bandpass integration on those templates, the output map is already smoothed. For a model like GNILC dust, using N_side 2048 on Jupyter@NERSC, it would be 40s to apply alm2map to the template, 10s and 10s for spectral index and dust temperature.

File sizes

Size for each component in single precision, multiply by 3 for IQU, multiply by 2 for double precision.

Nside	map (float32)	Alm (complex64) 2 Nside	Alm (complex64) 3 Nside - 1
32	49.2 kB	34.3 kB	74.5 kB
64	196.6 kB	134.2 kB	296.4 kB
128	786.4 kB	530.4 kB	1.2 MB
256	3.1 MB	2.1 MB	4.7 MB
512	12.6 MB	8.4 MB	18.9 MB
1024	50.3 MB	33.6 MB	75.5 MB
2048	201.3 MB	134.3 MB	302.0 MB
4096	805.3 MB	537.1 MB	1.2 GB
8192	3.2 GB	2.1 GB	4.8 GB
16384	12.9 GB	8.6 GB	19.3 GB